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PREFACE 



/'The main objective of this study is to develop a numerical model for 
routing water and sediment on small agricultural catchments. Part 2 of 
this report presents a detailed description of the models—The model/is 
developed on a general basis so that it may be applied to any agricultural 
catchment, and it can be used to simulate the effect of different land uses 
on the water and sediment yields from the modeled catchment. ^Paxt 

C h f * < 

presents results from validations of the model using several data sets 
including data from natural catchments. Input data and computer coding 
details are given in the Addendum to this report. 
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INTRODUCTION 


Alluvial streams are dynamic systems that continuously change their 
configuration and state in response to either changes in the natural 
environment, or perturbations introduced by man's activity. Frequently, 
these changes conduce to alteration of the stream-channel stability, which 
often results in channel migration and shoaling. 

Among the leading causes of channel instability are several that are 
intimately associated with land-management and conservation practices 
carried out on the upland areas. They are (a) clearing of land that 
removes the soil-protective and flow-retardant ground cover, which in turn 
leads to increased erosion and flood peaks; (b) installation of reservoirs 
for flood protection and irrigation control, which upset the water-sediment 
equilibrium downstream of those structures; and (c) excessive soil erosion 
resulting from uncontrolled sources. The combined effect is an aggregate 
flow of water and sediment coming from a variety of point and non-point 
sources within the upstream catchments. This aggregate yield acts as a 
time and space dependent loading on the streams draining the catchments. 
If this loading becomes quite different from that which the streams have 
adjusted to, the result is a breakdown in the stability of the channel 
system. 

The catchments contributing to the loading of any given channel system 
exhibit in general a great variety of soils, vegetation, and land uses. In 
order to effectively assess the impact of these catchments on the loading 
of the channel system, it is necessary to develop improved methods for 
predicting the effects of alternative land managements of those catchments. 
There is, therefore, a need for the development of mathematical models so 
that the hydrology of a catchment can be simulated and the effects of 
various management practices understood and predicted. In response to this 
need, the goal of the present study is the development of a prediction 
model for estimating sediment yield from agricultural catchments. In 
developing the model, the following specific objectives were considered: 
(i) estimate the amount of soil loss from specified soil-source units with 
homogeneous characteristics; (ii) estimate the amount of water and sediment 
transported out of the catchments through the principal drainage networks; 
and (iii) estimate the rate of channel aggradation and degradation along 
the flow system. The model is oriented towards the needs of the Corps of 
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Engineers for better means of assessing the impact of land-management 
practices on stream channel behavior. 

Many hundreds of papers have been written concerning studies on 
various aspects of hydrology. For this reason it is quite impossible to 
summarize the previous work that has led to the current understanding of 
the hydrologic cycle. Reference can be made to some of the existing 
comprehensive hydrology books (i.e., Chow, 1964), and to a number of 
American Society of Agronomy Monographs (Luthin, 1957; Van Shilfgaarde, 
1974; Hagan, Haise and Edminster, 1967) and reports (Pierre et al., 1966; 
Neilsen, Jackson, Cary, and Evans, 1972) that discuss current knowledge of 
soil-water-plant system. Excellent reviews of the progress made during 
recent years in several hydrologic subjects have been presented by Schaake 
(1975), Amerman et al. (1975), Johnson and Meyer (1975), and Nordin (1975). 
Only because of this accumulated knowledge is the proposed project even 
feasible. As a better understanding of the hydrologic cycle and the basic 
physical laws governing it have evolved, they have been synthesized into 
more rational and physically based models. Finally, the accessibility to 
high-speed digital computers has made possible the development of detailed 
comprehensive hydrologic models. 

Because the physical processes governing catchment behavior are very 
complicated, many past studies have utilized regression models. However, 
it is difficult to predict the response of a catchment to different land- 
management activities using regressional methods, because these methods are 
based on the assumption of time and space invariability. This assumption 
almost always fails to be valid in the case of natural catchments. 

A second type of models includes lumped parametric simulation methods, 
such as the TVA Continuous Daily-Streamflow Model (TVA, 1972). These 
models simulate the response of a given catchment by adjusting a number of 
coefficients, with little physical significance, using data collected under 
certain environmental conditions. The impossibility of relating those 
coefficients to a different set of environmental conditions, seriously 
restrict the use of these models for predicting the response of ungaged 
catchments. 

A different class of models embodies the distributed process 
simulation methods. These techniques use mathematical descriptions of the 
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basic hydrologic processes being modeled, and their interaction. In 
addition, this approach tends to minimize the number of adjustable 
parameters and, whenever possible, relates them to physical quantities that 
can be readily measured in the field. 

The Stanford model (Crawford and Linsley, 1966) was one of the first 
general models developed to simulate runoff from a catchment. It is 
basically a lumped-parameter model, although large, heterogeneous 
catchments can be subdivided into subcatchments if sufficient data are 
available to define model parameters. The model has gained widespread use 
and as a result has undergone numerous modifications. Holtan and Lopez 
(1970) have described the USDAHL-70 model of catchment hydrology. Although 
this model is basically lumped, a heterogeneous catchment can be broken 
down into smaller homogeneous areas. An attempt is made to incorporate 
spatial variability by dividing the catchment into land capabilities 
classes that correspond to uplands, hillslopes, and bottom lands. Dawdy et 
al. (1972) reported on a lumped-system model similar to the Stanford model 
which describes surface runoff from small catchments. TVA (1972) recently 
described a lumped daily-streamflow model with sixteen parameters, five of 
which require optimization. This model has been reasonably successful in 
predicting daily streamflows. 

A continuous distributed model is not yet available. However, several 
single-event distributed models that include part of the hydrologic cycle 
have been introduced since the pioneering works of Wooding (1966) and 
Woolhiser and Liggett (1967). Since then, a cascade of various sizes and 
slopes (Brakensiek, 1967; Kibbler and Woolhiser, 1970) or converging 
inverted cone-shaped surfaces (Woolhiser, 1969) have been used for 
geometric representation of complex topographies. The works of these and 
other investigators have let to the acceptance of the kinematic-wave 
approximation as an adequate model of shallow overland flow and flow in 
channels. The reductionist approach to watershed simulation was introduced 
by Huggins and Monke (1970). They employed a square grid for decomposing a 
complex catchment into elemental surface units. Most physically-based 
overland flow models used simplified lumped-system infiltration models. 
Smith and Woolhiser (1971) were the first to introduce a distributed 
infiltration model, derived from soil moisture flow theory, to calculate 
point infiltration rate, and therefore rainfall excess rate. The foregoing 









concepts have been incorporated, in one way or another, into more detailed 
models recently reported by Simons et al. (1975) and Smith (1976). These 
models use flow routing techniques based on the kinematic-wave 
approximation of the flow governing equations. Since its formulation by 
Lighthill and Whitham (1955), the kinematic wave approximation has received 
extensive application to catchment runoff modeling. This approximation is 
restricted by the assumption that the friction slope equals the stream bed 
slope, but it has been found to be applicable in many stream flows and in 
most overland flow situations. In addition, the kinematic-wave formulation 
admits an analytical solution by the method of characteristics (Eagleson, 
1970; Kibler and Woolhiser, 1970; Li et al., 1975b; Singh, 1975). This 
analytical solution has two main advantages over other numerical solutions. 
It eliminates the wave-celerity-damping and phase lag usually induced by 
numerical schemes; and, in addition, results in faster computational 
procedures. In spite of these advantages, applications of this analytical 
solution have been restricted in the past to catchment models with a high 
degree of geometric abstraction. The reason is the formation of kinematic 
shock waves (Lighthill and Whitham, 1955; Kibler and Woolhiser, 1970; 
Harley et al., 1970; Whitham, 1974). Formally, innumerable shock waves can 
be generated during the routing process, as a result of the time and 
spatial discretization of precipitation and the physical characteristics of 
the catchment. In the past, the existence of these shocks has frequently 
been ignored by using approximate numerical techniques. This practice, 
however, may not be considered as valid particularly when the foundation 
and the physical relevance of the kinematic wave approximation is under 
investigation. It is well known that shock formation is intrinsic to the 
hyperbolic equation governing kinematic theory. Further, they are 
considered to be the manifestations of higher order effects such as 
formation of monoclinal flood waves, bores, etc. These discontinuities 
play important roles in the dynamics of hydraulic systems and an ad hoc 
smoothing by numerical means does not necessarily make the theory look 
better. The model described in the present report introduces a new 
solution to the kinematic approximation, which retains the dynamic effects 
of the shocks by routing the discontinuities as they appear. Certain 
simplifying assumptions are made which permit closed form solutions and an 
efficient numerical algoritlim, based on the method of characteristics. The 
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resulting procedure, called an approximate shock fitting scheme, preserves 
the effect of the shocks without the usual computational complications and 
compares favorably with existing finite difference solutions (Borah et. 
al., 1980). 

Different types of sediment production models have appeared widely 
dispersed in the technical literature. Reference can be made to a recent 
review presented by Heinemann and Piest (1975) and to a publication of the 
Agricultural Research Service (1975). Several regression equations for 
predicting gross soil erosion have been proposed. The most commonly used 
among these is the so-called universal soil loss equation (USLE) proposed 
by Wischmeier and Smith (1965). Other equations of similar nature have 
been developed by Musgrave (2947), and Gottschalk and Brune (1950). In 
these equations, the soil loss rate is correlated with storm, land, and 
vegetation characteristics. Such equations are applicable on seasonal 
basis or longer. Also, they do not take advantage of the physical 
processes occurring within the catchments; hence, it is not possible to use 
them on large, complex basins. Williams (1972) modified the USLE to make 
it applicable for predicting storm sediment yields. Onstad and Foster 
(1975) combined a different modification of the USLE with the USDAHL-70 
catchment model to predict sediment yield for single storms. They applied 
their model to two small catchments with limited success. 

The first physically-based sediment yield model was reported by Negev 
(1967). This model uses the Stanford model for the water phase, and takes 
into consideration rainfall soil splash, entrainment by overland flow, and 
rilling and gullying, along with separate channel transport of fine and 
coarse sediment. Sediment production is evaluated in terms of power 
functions of water discharge containing a number of parameters that must be 
calibrated. A modified version of Negev's model has been incorporated in 
the Agricultural Runoff Management (ARM) model recently reported by 
Donigian and Crawford (1976). The aforementioned models of Simons et al. 
(1975) and Smith (1976) also incorporate the capability of describing 
sediment movement on a catchment as a time and space distributed process. 
The structures of these two models are similar; however, there are 
differences in numerical techniques and functional relationships. The 
sediment movement is described by linking the excess-rainfall flow 
equations to the sediment continuity equation, with relations describing 
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sediment detachment and transport capacity at any point on the surface or 
in a channel. A similar structure has been incorporated in the erosion and 
sedimentation component of the present model. In addition, sediment is 
routed using a sediment characteristic scheme that takes advantage of the 
efficient analytical solution mentioned above. 

Part 2 of this report provides a detailed description of the model. 
Part 3 discusses the validation of the model on sets of laboratory and 
field data. Input data and coding details are given in Addendum 1. This 
report is based on material presented in an earlier study by Borah (1979). 
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MODEL FORMULATION 


The model basically consists of two intertwined models: one 
describing the hydrology of the basin; the other describing the associated 
erosion and sedimentation processes. It simulates the movement of water 
and sediment as a time and space distributed process, and it has the 
ability of distinguishing between overland and channel flows. The 
catchment is regarded as consisting of a mosaic of individual subcatchments 
interconnected by channel reaches. The model can thus be regarded as a 
cascading process in which the output of one or more subcatchments becomes 
the input to another subcatchment or channel reach. In mimicking the 
overland movement of water and sediment, the model simulates processes of 
interception, infiltration, runoff, detachment, transport, and deposition 
of sediment. The water and sediment reaching the streams are then routed 
through the channel system, and the rates of channel aggradation and 
degradation are computed. The basic structure of the model is graphically 
illustrated in Fig. 1. The details of the model components are given in 
the following sections. The applicability of the model is restricted to 
catchments where the streamflows are ephemeral, the subsurface flow and 
ground water movement are not significant, and the kinematic wave 
approximation for flow routing is valid. The computer program has been 
written assuming a uniform distribution of rainfall. However, the program 
can be easily modified to accomodate any other aerial rainfall 
distribution. 

The catchment is segmented into subcatchment and channel reaches to 
account for the lack of uniformity in terrain, soil, and land use 
characteristics in most natural catchments. Within each of the segments 
these characteristics are treated as being uniform. The subcatchments are 
replaced by sloping rectangular areas with representative length, slope, 
width, soil, and vegetative characteristics. The channel segments are 
described by representative cross-sectional shape, slope, length, and 
roughness. Gravity flow logic is used to determine the computational 
sequence as explained in Appendix I. The input data required by the model 
includes storm characteristics, geometry data, vegetative cover data, soil 
data, and water and sediment routing data. Details on input data 
preparation are given in Addendum 1. 
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Fig. 1. Structure of the model 
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2.1 INTERCEPTION. NET RAINFALL RATE 

The rainfall excess is the rainfall contributing to the water flowing 
over the surface of an overland unit. This is the resultant rainfall after 
incurring the losses due to interception, evaporation, transpiration and 
infiltration. At the beginning of a rainfall event, some rainfall is lost 
due to interception at the canopy and ground covers and thus interception 
is the first concern in computing rainfall excess. The interception 
component adopted in this model is based upon the approach proposed by 
Simons, et al. (1975). 

The canopy cover and the ground cover are the two major features which 
influence the motion of raindrops before reaching the ground surface. The 
canopy cover and the ground cover are represented by the canopy cover 
density and the ground cover density, respectively. The canopy cover 
density, D c , is defined as the ratio of the area covered by trees to the 
total area, and the ground cover density, D^, is the ratio of the ground 
area covered with litter, rock, grass, etc., to the total area. The canopy 
cover and the ground cover densities are assumed uniform on a flow unit. 

Net rainfall is the quantity of rainfall reaching the ground after 
passing through those covers. Without any obstruction like canopy cover or 
ground cover, the rain reaches the ground surface at the same intensity 
without any recognizable loss. In the presence of trees, a portion of the 
rainfall is stored in the canopy and the remainder passes through the trees 
as throughfall and stemflow (Zinke, 1965). In the presence of ground 
cover, a portion of the throughfall is intercepted and the remainder 
finally reaches the ground as net rainfall. 

Let I be the rainfall rate (or intensity) at time t at a level above 
the tree canopy as shown in Fig. 2. Let I be the rate at which rain is 
being stored in the canopy at time t. Then, the rainfall rate under the 
tree canopy is reduced to the throughfall rate and the weighted average 
throughfall rate is 

I = I - D I , (1) 
o c c 

where stemflow has been neglected. Similarly, let I be the rate at which 
rain is being stored in the ground cover at time t. Then the weighted 
average net rainfall rate reaching the ground is 

I = I - D I (2) 
n o g g 
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2 Control volumes for tree canopy and ground cover 
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According to Horton (1919), the total interception equals leaf storage 
capacity plus evaporation loss during the storm. Zinke (1965) indicated 
that . . usually for a storm, there is an initial period during which 
the vegetation cover is wetted and a so-called interception storage 
capacity is satisfied. This is followed by loss from this storage, and the 
loss is dependent upon the evaporation opportunity during the remainder of 
the storm." 

Let be the interception storage capacity of a tree canopy per unit 
area of tree canopy, the interception storage capacity of the ground 
cover per unit area of ground cover, E the mean evaporation rate form the 
interception storages, S c and the ratios of the evaporating surface area 
to the horizontal projected area for a tree canopy and for a typical ground 
cover, respectively, and I the initial interception storage content which 
is defined as the ratio of the initial storage to the interception storage 
capacity. Accordingly, 


I c = I, if , / K*)d t S (1 - I s ) V c , 


I = ES , if , f I(t)dt > (1 - I ) V 
c c J v s c 

o 


Similarly, 


I = I , if , J I (T.)d < (1 - I ) V 
go’ to s' g 


I = E S , if , J I (t)dt > (1 - I ) V 
g g o ° S 8 
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For easy handling of data, a parameter r v is introduced, defined as 

the ration of the interception storage capacity of a typical canopy cover 

to that of the ground cover. Therefore, 

V = r V 
c v g 

and 

S = r S 
c v g 

Eqns. 3 and 4 can be rewritten in the discrete form as 

I = I , if I IAt < (1-1 )r V , (7) 

c S V g ’ ' ' 

I = E R S , if I IAt > (1-1 )r V . (8) 

c v g s' v g ' 

Similarly, one obtains 

I = I , if I I At < (1-1 )V , (7’) 

go’ o s' g 

I = E S , if I I At > (1-1 )V . (8’) 

8 g o s' g 

Eqns. 1, 2, 7, 8, V and 8* are used to compute the discretized net 

rainfall rate. 

2.2 INFILTRATION. RAINFALL EXCESS RATE 

Water reaching the ground surface at the beginning of a storm passes 
through into the soil because in nearly all cases the initial infiltration 
capacity rate is greater than the initial net rainfall rate. In that case, 
the infiltration rate is equal to the net rainfall rate and there is no 
runoff. Once the net rainfall rate exceeds the infiltration rate, the 
excess rainfall either runs off or is stored on the ground surface as 
depression storage and thus ponding takes place. The time elapsed until 
the beginning of excess rainfall is known as the ponding time. 
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It is very difficult to describe mathematically the process of runoff 
originating from rainfall excess, because, very little is known concerning 
the magnitude of depression storage. Meaningful observation of depression 
storages are not easily obtained. Thus, depression storage is usually 
combined with interception and treated as an initial loss with respect to 
storm runoff (Linsley et al. 1958). For simplicity, the depression storage 
is omitted in this model, but implicitly is included in the interception 
storage capacity described in the previous section. 

Thus, the water balance equation at the ground surface is 


I = I - f. 
e n 1 


(9) 


in which I is the rainfall excess rate, I is the net rainfall rate, and 
e * n ’ 

f. is the infiltration rate, 
l 

The ponding time and the infiltration rate are computed from the 
infiltration model developed by Smith and Parlange (1978). The model was 
developed from a simplified solution of the equation for one dimensional 
diffusion of water under gravity, by imposing the rainfall pattern as 
flux-type boundary condition. The equation and the general initial and 
boundary conditions are 


90 _ 9_ 
at “ 9z 


<D §i ) 



( 10 ) 


and, 


t = 0, z > 0, 0 = 0^ ; (11a) 

0 < t < t p , z = 0, K-D~ = I n (t) , (lib) 


where 0 is the soil water content; z is the vertical distance from the soil 
surface; D is the moisture diffusivity; K is the hydraulic conductivity; t p 
is the ponding time; and 1 n (t.) is time-varying net rainfall intensity 
pattern. Eq. (10) is also generally written as 


9z 9_ ,_90. _ dK 
at 90 1 9z' d0 


( 12 ) 


by using an elementary identity of differentials. The primary assumption 
used in the derivations is that D varies rapidly with 0, which implies that 
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water flux within the soil varies little with relative position. 
Integrating Eq. 12 with the condition 11b, and using an expression of 
conservation of mass, Smith and Parlange (1978) obtained the solution 


t 8 

; I dt = S (0-0;) 


I -K(0) 
n 


where the subscript 1 refers to conditions at the surface, and i refers to 
initial conditions. When ponding occurs 8^ = 8 g (saturated water content), 
and the left-integral upper limit becomes t = t . 

By further assuming that K varies slowly in the region near 8 s> or if 
K is much larger than I , Eq. 13 may be approximated by 


S v rdt = B(0.) (I - K ) 
Q i np s 


where I is the net rainfall rate at ponding time, and 6(8^) is a 
parameter dependent only on soil type and 8^. B is theoretically 
identified as a function of sorptivity and it can be roughly estimated from 
B = S 2 /2, where S is the soil sorptivity. Eq. 14 is used in estimating the 
ponding time, t . Where I n (t) is a continuously differentiable expression, 
an analytical expression for t results from Eq. 14. Engineering use of 
this expression is quite simple. As a rain pattern progresses, the value 
of iv* is calculated and compared with the value of the right hand side 
of Eq. 14, with I = 1^. This expression is an inequality only up to the 
time at which ponding occurs. 

For the same conditions on K given above, the expression for the decay 
of infiltration capacity for t > t was derived as was the expression for 
ponding time. A more general expression for surface conditions related to 
Eq. (13), which reduced to Eq. 14, is 


o (8-8.)K 
t . n 


df = J (0-0.)KdY S I dt 
4* 1 o n 


Here, the variable of integration is f, rather than 8. This expression 
applies for all times. The value of 4* at the surface is taken equal to 
zero for t > t , although a surface depth could be taken into account if 
desired, the effect of which is, in fact, small. 


• 'V* 










Differentiating Eq. 15 with respect to t yields 


- r IT ° s (e " e i ) TT^kT 5 ’ ™ = ? (Q-ejKdf 

n dt f. 1 U n KJ 4<. 

x i 


(16) 


Taking t > t^, the surface flux is f(t) < I n (t), and integrating from t to 
t yields 

° ( 0 - 0 ^ I n (f-K) y v 

A l f nr^ *" 'nf-Ttr 1 = 

•F np np-K 


(t-t ) J (e-e.)Kd'f 

P y L 

i 


(17) 


Here A depends slightly on the rainfall rate pattern. From Eq. 14 


A = B = (I -K ) y 1 dt 
np s J n 
r o 


( 18 ) 


The function of K in the left integrand of Eq. 17 can be replaced by its 
value at saturation, K g , and combined with Eq. 18 to eliminate A, yielding 


t 

K s (t-t ) = J P rdt {~ 

' " o s 


I -f I -K (I -K )f 

_Q£_s o n r_5£__§—1} 

f-K K l (f-K )I u 

s s np 


(19) 


Eq. 19 is solved numerically to compute the infiltration decay, f(t) 
for t > t . In this model, Newton's method is applied and according to 
this method, the infiltration rate at the ith time interval and kth 


iteration step is 


= f 


F < f i,k-1> 


i.k i i,k-i F'd. ^) • 


( 20 ) 


The function F(f) and its derivative, F'(f), are 


I -f I -K 


I -K 


F(f) = K (t-t ) - J* I dt {JE- f-i tn (-f-J 7^)1 = 0, (21) 

0 s s np s 
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and 



To start the computation, the following initial values are assumed: 


f 

0,0 


f(tp) = I 


np 


(23a) 


f . 
1,0 


(23b) 


Due to the lack of infiltrometer data in most of the catchments, the 

parameter B in Eq. 14 is presently estimated from the results presented by 

Smith and Parlange (1978), Fig. 3. For simplicity, the results for Colby 

S. L. (swelling) were used in the tests discussed in this report. With the 

known information of initial soil moisture deficit, (0 -0.) and the 

si 

saturated hydraulic conductivity K g (cm/sec), the value of B = S 2 /2 

(cm 2 /sec) is estimated from the quasi-linear functions plotted in Fig. 2. 


2.3 WATER ROUTING 

The motion of water moving either as overland flow or channel flow, is 
governed by the equations of mass continuity, momentum balance, and flow 
resistance. The flow routing scheme described in this report is based on 
the kinematic-wave approximation to the general equations of motion. The 
following sections present the basic assumptions underlying the kinematic- 
wave theory, along with the essential points of the analytical solution. 
This solution is then used to investigate the conditions under which 
kinematic shock waves may be expected. Subsequently, procedures for proper 
shock fitting are discussed. Simplifying approximations are made which 
allow closed form solutions and preserve the effects of the shocks. The 
resulting approximate shock fitting scheme is compared with an existing 
implicit finite difference solution. The accuracy and efficiency of the 
new scheme are illustrated in Part 3 by computing a variety of unsteady 
flows, ranging from simple cascades to complex natural catchments. 
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9 - 9 . INITIAL SOIL MOISTURE DEFICIT 
s ® 


Fig. 3 Quasi-linear relation of the second power of soil sorptivity to 
initial moisture deficit (after Smith and Parlange, 1978) 



2.3.1 Kinematic Wave Approximation 

This approximation is based on a simplification of the Saint-Venant or 
shallow-water equations governing unsteady free-surface flows. The 
underlying assumption is that energy gradients due to local and convective 
accelerations are negligible in comparison with gravitational and 
frictional effects. The momentum equation thus become: 


S 0^ 


(24) 


where and S^. are the bed slope and friction slope, respectively. The 
continuity equation for water is as usual: 


9a + ag = 

3t 9x 




(25) 


in which A is the flow cross sectional area, Q is the flow rate of 
discharge, x is the downslope position, t is time, and q^ is the rate of 
lateral inflow or outflow per unit length of stream. 

Any suitable law of flow resistance can be used to express Eq. 24 as a 
parametric function of the stream hydraulic parameters. A widely used 
expression is as follows: 


Q = aA n 


(26) 


where or and n are parameters related to channel (or plane) roughness and 
geometry. Obviously, a and n are functions of time and position. The 
space dependence can be removed by making the stream hydraulic properties 
and the lateral inflow piece-wise uniform in space (Fig. 4). The flow 
region is segmented into a network of different elements with properties 
remaining constant within each element, but varying from element to 
element. This approach, known as a kinematic cascade, was introduced by 
Brakensiek (1967b) and latter elaborated upon by Kibler and Woolhiser 
(1970). Similarly, the time dependence of those coefficients and the 
lateral discharge can be eliminated by assuming each piece-wise constant in 
time (Fig. S). Then Eq. 23 and Eq. 2b can be solved for each cascade 
element subject to given initial and upstream boundary conditions (i.e., 
upstream inflow rate, Qy)- The subscripts U and D will be used to indicate 
upstream and downstream conditions, respective 1 y. 
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2.3.2 Analytical Solution 

Eqns. 25 and 26 can be solved using the method of characteristics 

(Eagleson, 1970). The basic steps of this method are given below for 
completeness. Assuming that a and n are constant on a given cascade 
element k and combining Eqns. 25 and 26 yields 


9A 

at 


+ an A 


n-1 3A 
9x 




(27) 


The total differential of A yields 


8A 

at 


dt + dx = 
dx 


dA 


(28) 


Since the derivatives in Eqns. 25, 27, 28 do not exist along the 
characteristic paths, the determinant of the coefficient matrix 
corresponding to those equations must vanish at points on the 
characteristics. This requirement thus yields 


dx 

dt 


= omA 


n-1 


1 (1 - -) 
na 11 Q n 


(29a) 


or, 


dt 

dx 


anA 


n-1 


Obtaining the invariants of the solution gives 


dA _ 
dt q £ 


(29b) 


(30a) 


or, 


dA 

dx 


anA 


n-1 


(30b) 


The sets of Eqns. 29A, 30A or 29b, 30b represent the characteristic 
form of the solution to the kinematic-wave approximation. They show that 
kinematic waves possess only one system of characteristics. Accordingly, 
kinematic-wave routing cannot be used in situations where there are 
downstream flow controls. Integrating Eqns. 30a and 30b with the initial 
and boundary conditions 
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Aj = A(x,ty), t Q = 0, 

(31) 

Ay — A(Xy,t), x Q - 0, 

(32) 


results in the following respective expressions for the flow conditions 
along any characteristic t = t(x) (Harley et al., 1970): 


and 


A(x,t) = A + f q.{x, £(x)] d£ , 
l 0 


A(x,t) = J A o + h f do) n , 


(33a) 


(33b) 


or 


Q(x»t) = aA[j + J q £ (n,t(n)l dr). 

x o 


(33c) 


Substituting Eqns. 33a and 33b into Eqns. 29a and 29b, respectively, and on 
integration yield 


t t' 

x-x Q = an / {/ q £ [x,£(x)J d£ + A Q } 

C 0 t 0 


n-1 


dt', 


(34a) 


and 


1-n 


t_t 0 = k f { a f q £ ln ’t(n)J dn + aJ} n dx’ . 
x o x 0 

In these equations, A Q represents the initial flow area, Aj 


(34b) 


along the 

characteristic (shown in Fig. 5), and the upstream inflow area, Ay = 
1 1/n 


[Q u (t)/o] 


along characteristics like Cj, C £ , and C^. The above 


integrals are functionally integrable and yield simple expressions when 
q £ (x,t) is either an explicit function of x and independent of time, or 
vice versa. In most applications, only discrete distributions of lateral 
inflows are available. For instance, rainfall intensities and/or 
infiltration rates are usually regarded as lateral inflows in hydrologic 
models. In these cases, it has been a common practice to use lateral 
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inflows having a piece-wise uniform distribution in space (Fig. 4) and a 
piece-wise constant variation in time (Fig. 5). This form of lateral 
inflow distribution is also adopted in this paper. Thus, assuming q^(x,t) 
remains constant within small space and time increments the above 
expressions may be explicitly integrated. Discretizing Eqns. 33 and 34 
between two points and (x^tj) on a characteristic path (Fig. 

5) yields 


A. . = A. _ . , + q n . At. , 

i.J i-i»j-i £,J J 

1 

A. . = (A. n . . + Ax. ) n , 

i.J i-l,J-l a i 


Q i,j = Q i-l,j-l + **i ’ 


and 


Ax. = 
i 


[ (q„ . At. + A. . ) n - A. n . , J , 

q„ . £,J J x-1,j-l i-l, j-l 

*»J 


At. = -1— [(^J- Ax. + A, n . _)" - A. , . J , 
j a i x-i, j-r ’ 


1 

,n 


(33'a) 
(33'b) 
(33'c) 
(34'a) 

(34'b) 


where Ax = x. - x. ,, At. = t. - t. , and q» . is the rate of lateral 

inflow assumed constant over the time step At. and the space increment Ax.. 

J 1 
The last two equations are not defined when q„ . =0; in this case Eqns. 

“, J 

29a and 29b yield 


a »n-l 

Ax. = an A. , . , At. 
i i-l,j-l j 


and 


A . 1 ,1-n . 

At. = — A. . . . Ax. 
j an i-l,j-l i 


134’c) 


(34'd) 


Eqns. 34' are used to trace the characteristic path by considering 

either Ax^ or At^ as a dependent variable and choosing a suitable value for 

the other. This increment must be chosen so that the lateral inflow can be 

assumed constant within the intervals At. and Ax.. Eqns. 33' are used in 

J i 

turn to compute the flow conditions on the same characteristic. Eqns. 33' 
and 34' are explicit and independent of each other, and therefore, they may 
be used in several different combinations. In this paper, the following 
scheme is used for computational efficiency. Since the lateral inflow is 
uniform along a cascade element and piece-wise constant with time, the best 
choice is to take the time increment At. as the independent variable. This 


1.28 








time-increment is used in Eq. 33'a to compute A^ ^, and Eq. 26 is used to 

obtain Q. .. Then Eq. 33'c is solved for Ax.. These steps constitute the 
1 > J i 

scheme used to trace a characteristic across the cascade element and to 

compute the flow conditions along its path. 

In general, a characteristic does not intersect the downstream 
boundary exactly at the beginning or end of a time step. A modification of 
the proposed scheme is used to determine the time of intersection. The 
last space increment Ax^ = x^_j, where x^ is the length of the cascade 
element, is substituted in Eq. 33'c to compute and then Eq. 26 is solved 
for Ap. This value is then introduced into Eq. 33'a to solve for the time 
increment, At.., and the time of intersection t_ = t„ , + At VI . The 

discharges existing on all the characteristic paths at the time of their 
arrival to the downstream boundary define a discrete outflow distribution. 
This discrete distribution is smooted out by interpolation to obtain a 
continuous outflow hydrograph. The problem of hydrograph computation for 
any given cascade element is, thus, completely specified once the initial 
flow areas along the element at time zero (Eq. 31) and the inflow 
hydrograph, coming from the upstream element (Eq. 32), are known. 


2.3.3 Formation of Shock Waves 

The foregoing analytical solution remains valid as long as the 
characteristic paths do not intersect each other. When this occurs, the 
preceding theory does not give a unique solution since there is more than 
one characteristic passing through the intersection point, causing flow 
discontinuity. These discontinuities have properties analogous to those of 
shocks that arise in gas dynamics theory and, by analogy, are generally 
known as kinematic shock waves. Eq. 29a shows that the celerity of a wave 
moving along a characteristic path is a function of the flow area A. This 
dependence on A produces a nonlinear distortion of the wave as it 
propagates. Waves with higher values of A travel faster and finally 
overtake lower ones leading to the breaking of the wave profile as 
illustrated in Fig. 6. In this figure the lateral inflow has been ignored 
for simplicity and, thus, the flow areas do not change along the 
characteristic paths (Eq. 30). The wave profile at any instant t = tj > 0 
is obtained by projecting the corresponding values of A on the verticals 
passing through the intersections of the characteristic paths and the line 
t = t,. The wave breaking begins at the time t = t_ when the free surface 
profile first develops an infinite slope. In the x-t plane this breaking 
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coincides with the intersection of two characteristics, where a single¬ 
valued flow area no longer exists. This sudden change in flow area is 
interpreted as the initiation of a shock. After its formation, the shock 
proceeds downstream with its own velocity, as discussed in the next 
section. 

However, a shock will be formed only if, given any two consecutive 
characteristics, the first one is less steep than the following one. For 
example, Fig. 5 shows the characteristics and passing through the 
points E and D at the same time t'. Then, the necessary condition for 
shock formation is that 


dx (1) 
f—1 
dt E 


[ —] 
l dt J 


( 2 ) 

D 


(31) 


where the superscripts refer to the consecutive characteristics and Ca¬ 
using Eqns. 26, 29a and 33a, Eq. 35 can be expressed as 


1 


1 


1 


r _ ( 2 ) , n k f . (1), n k . , "k A . 

lQ U,k 1 ‘ lQ U,k 1 ( V q 2,k At j + 1 


(36) 


in which the subscript k indicates flow variables associated with the kth 

cascade element and At.^, = t.^, - t.. This inequality indicates that 

J + 1 J + 1 J 

characteristics emanating from the line t = 0 cannot intersect for either 
dry or uniform initial conditions. It is also clear from Eq. 36 that 
characteristics originating during steady upstream inflow or along the 
falling limb of the upstream hydrograph (i.e., characteristics in Fig. 
5) will not intersect, unless the flow law changes from turbulent to 
laminar. These two zones are then free from shock formation. This does 
not mean that shocks originating upstream may not propagate into these 
zones. On the other hand, the characteristics originating along the rising 
limb of the upstream hydrograph may intersect depending upon the relative 
magnitude of the variables appearing in Eq. 36. For this reason, the 
domain bounded by the characteristics emanating from boMi ends of the 
rising limb can be regarded as the probable shock-forming zone. The 
solution in the shock-forming zone is influenced by the upstream boundary 
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conditions of the kth-element which, in turn, depend on the initial 
conditions on the element k-1. For example, applying Eq. 36 to the 
characteristics passing through the points M and N indicated in Fig. 5, and 
using Eqns. 29a and 33c, one obtains 


k-1 

(JLJ.) 

v /M ' 


1. 

n, 


k-1 




’t.k'Vi 


(37) 


This inequality shows that for uniform topography (a, = a , n, = n ) a 

K“1 K K"i K 

shock will be formed in element k whenever q^ k ^ >q^ whereas, under 
conditions of uniformly distributed lateral inflow (q 0 = q ), a shock 

will occur wherever a kl > a k for n kl = n k and n k J > n k for a k = a k * 

Harley et al. (1970)., Kibler and Woolhiser (1970), and Li et al. 

(1976a,b) also studied the foregoing necessary conditions. The latter also 
discussed the sufficient conditions for two given characteristics to 
intersect within the boundaries of the cascade element. They suggested 
taking two consecutive characteristics far apart enough so as to avoid 
their intersection within the element (i.e., characteristics C' and C" in 
Fig. 5). However, once the necessary condition is met, the shock wave S 

may form within the cascade element and travel all the way to the 

downstream end (this is discussed in the next section). Therefore, not 
considering the solution domain between C' and C" serves only to ignore the 
possible existence of the shock but does not really avoid it. Given that 
the shock is an intrinsic part of the solution to the kinematic-wave 
approximation, a better routing procedure is one that considers the 
presence of the shock and its effect cn the outflow hydrograph. Such an 
approach is discussed in the next section. 


2.3.4 Approximate Solution With Shock Fitting 

After breaking occurs the kinematic wave theory ceases to be a 
strictly valid description of the physical process, because the flow area 
is inherently single-valued. However, the foregoing formulation can still 
be utilized by allowing discontinuities in the solution. Lighthill and 
Whitham (1955) proposed replacing the shock, as a first approximation, by a 
discontinuous wave (BC, Fig. 6) that produces the appropriate variations in 
discharge and flow area as it moves downstream. These authors derived the 
velocity of the shock by considering the continuity of flow and the rate of 
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discharge crossing the shock front. Applying their result to an arbitrary 
points (x,t) on the shock path resulted in 


u(x t) = dx Q Jx 2 t) 

*• ’ J dt ,b 


Q (x,t) 


(38) 


A (x,t) - A (x,t) 


where the superscripts b and a indicate flow conditions behind and ahead of 
the shock, respectively. While the shock is moving downstream, 
infinitesimal waves traveling along characteristic paths will join the 
shock from ahead and behind continuously modifying the shock strength and 
velocity. Therefore, the locus of all the points, where those 
characteristics intersect each other, defines the path of the associates 
shock. 

The important steps involved in routing a shock include finding the 
position where it originates and then tracing the shock downstream. An 
analytical method for finding the position where the shock originates was 
discussed by Whitham (1974). Because of its computational intricateness, 
his approach was not used in the present study, but, instead, the following 
simpler technique was developed. It is clear from Eq. 36 that once it is 
satisfied, the characteristics tend to converge. Thus, a shock may 
originate within the same element or in one of the following cascade 
elements, depending upon the magnitudes of a^, n^, and q^ ^ as they appear 
in Eq. .36. In this paper, the problem of finding the exact location of the 
shock origin is avoided by discretizing the upstream inflow hydrograph. By 
doing this, small artificial shocks are introduced all along the rising 
limb of the hydrograph, but they are routed across the characteristic plane 
only in those zones where Eq. 36 is satisfied. The remainder of the inflow 
hydrograph is routed according to the continuous solution from the 
characteristics. Thus, the shock fitting approximations introduced in this 
paper influence only those zones where shocks form, whereas the accuracy of 
the analytical solution is preserved in the other zones. The procedure is 
illustrated in Fig. 7, where individual characteristics are traced one at a 
time, starting at the upstream boundary of the characteristic plane. They 
are assumed to emanate from the half-time levels t^, tj + ^, t 2+V etc '* 
where the continuous and discretized upstream-hydrographs coincide. The 
shock-forming condition, Eq. 36, is tested at each new time step. The 


1.33 






















discharges and Qy^ appearing in this relationship are replaced by the 

discretized inflow rates at the new and previous time steps, respectively. 

In addition, the product ^At^ + ^ is replaced by the integral of the 

lateral inflow hydrograph over the interval t. , i t $ t. , , to account for 

J - "! J + "5 

the fact that the characteristic paths arise from t. , instead of t. 

J +! s J + l 

Hence, if j = 2, for instance, Eq. .36 is given by 


i_ J_ I_ 

(QU^) 0 " - (Qy,)^ > («/“ F 2 , 

where F 2> Qy 2 > an< * Qy 3 are defined in Fig. 7. When Eq. 36 is not 
satisfied, the corresponding characteristics are routed one at a time using 
Eqns. 33 and 34, as illustrated by C^, C^, and C in Fig. 7. On the other 
hand, if Eq. 36 is met, an artificial shock is introduced at the time level 
t^ and then routed using the scheme presented below. 

The problem of tracing the shocks in the characteristic plane has been 
discussed in detail by Whitham (1974). It consists essentially of 
determining the characteristics intersecting each other on the shock path, 
and the location of the point of intersection. The solution to this so 
called three-point problem, requires the simultaneous solution of the 
characteristic Eqns. 34 and the expression for the shock velocity, Eq. 38. 
This method was used by Kibler and Woolhiser (1970), who developed an 
iterative scheme to route shocks generated on a three-plane cascade under 
constant lateral inflow. In principle, this approach can be extended to 
include more complex geometries with varying lateral inflow. However, the 
computational details become excessively complicated, particularly when 
shock interactions occur. Therefore, an alternate approach is developed in 
this paper by restricting the above artificial shocks to small amplitudes 
(weak shocks). Consider the characteristic pair C^, intersecting at the 
point P on the shock path (Fig. 7), and a similar pair C^, joining the 
shock at an arbitrary point B. If the shock is weak, the flow areas 
carried by and (ahead of the shock) will not differ significantly 
from each ether. Similarly, and will have approximately equal flow 
areas behind the shock. Since B is arbitrary, and Clj may represent any 
of the infinite pairs of characteristics joining the shock (behind and 
ahead of it) along its trajectory. It is thus reasonable to postulate that 
the intersecting characteristics defining the shock path do not deviate 
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significantly from it. Then, as a further expedience, we assume that the 
shock path is defined by two sufficiently close (overlapping) 
characteristics, one representing the characteristics ahead of the shock 
and the other representing the characteristics behind the shock. These two 
characteristics will carry flow areas (depths) ahead and behind the shock 
as given by Eq. 33a or 33b. The common path of these characteristics, as 
well as the shock, is obtained by substituting Eqns. 33a and 26 or Eqns. 
33b and 33c into Eq. 38; after integration they yield 

t , t' 

X-X 0 = ~ET~a $ * lA 0 + q £ ( x »U x )) d£) n " IAq + 

V A 0 t 0 t 0 


/ q„( x >4(x)) d£] n } dt' , 


(39a) 


XX — 

t-tn = TT-^T / {I(A«) n + l S q 0 (n,t(n)) dp]" - [(A*)“ + 


x o 


l • / q £ (n,t(n)) dnl n ] dx’ , 
x o 


(39b) 


where (A^, Q q ) and (A^, Q^) are the known flow conditions ahead and behind 
the shock, respectively, at a given point (Xg,t^) on the shock trajectory. 

Consider now any two consecutive points B(x^ ^.t^^) and D(x^,t^) on a 
shock trajectory S (Fig. 7). Applying Eqns. 39 between these two points 
and discretizing the results, gives 


Ax. = 


(n+l)q. 


uv t r A h.i-i l " n 


'V'-j + 




At. = 


J (n+l)q £ (Q 


an i, i , , .,b ,i 

b--: u <r * ( A i-,,j-,) 

i-1, j-1 


a > /A a ,n, n , A b .n+1 

l T **i + (A i- lf j-i ) 1 - (A i-i, j -i ) + 


(A a ) n+1 } 


(40b) 

















in which the lateral inflow is taken as constant over the finite increments 


Ax. and At.. Eqns. 40 are not valid when q. = 0. In this case, a discrete 
i J * 

approximation of Eq. 38 is used. Similar to the characteristic wave 
routing, At^ is taken as the independent variable. Eq. 33'a is used to 

compute A 3 . and A b . and using these two in Eq. 26 0 a . and Q b . are 

^*.1 i i J i) J 

determined. Then Eq. 40a is used to compute Ax^. These steps constitute 

an explicit scheme that is used to route the shock like a characteristic 

wave. Similar to the characteristic wave, the scheme is modified near the 

downstream boundary. Here, Eqns. 33'c and 26 are used to compute the flow 

conditions ahead and behind the shock, and Eq. 40b is used to find the 

exact arrival time of the shock. Although this scheme is based on an 

approximation that deviates from the uniqueness of the solution, its use 

does not significantly affect the accuracy of the computations. In the 

foregoing discussion the shock fronts have been considered mathematically 

as flow discontinuities fitted into regions where multiple values of the 

solution exist. Physically they will not be discontinuities; instead, the 

fronts will have finite lengths induced by diffusive effects as well as by 

breaking of the waves. In many instances, the front will take the smooth 

shape of a monoclinal flood wave (Whitham, 1974). Because of these 

smoothing tendencies, and for the sake of computational expedience, the 

shock front can be approximated by a linear profile as shown in Fig. 7 

(segment EF). Thus, the flow condition at the shock itself can be computed 

as 

A(x,t) = | (A a (x,t) + A b (x,t)]. (41) 

This could be viewed as a method to make an approximately continuous 
solution from one which was made discontinuous by discretizing the inflow 
hydrograph. The preceding routing approach, called herein a propagating 
shock-fitting (PSF) scheme, permits the use of essentially the same 
numerical procedure to route both characteristic and shock waves. This 
scheme is particularly efficient when the initial flow conditions are 
either dry or uniform, and only the outflow hydrograph is of interest. In 
these cases the calculations are performed only in the x-t subdomain 
bounded by the lines C. and C (Fig. 7) where C is the characteristic 
reaching the downstream boundary at the end of the routing interval. All 
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the characteristics emanating from the line t = 0 are parallel to C^. 
Hence, the outflow conditions at times t^, t^, t^,... are equal to those 
computed at the points 1, 2, 3,... on C^. Within the above subdomain the 
step-by-step integrations along the wave trajectories use simple algebraic 
relationships. Moreover, values computed at interior points do not have to 
be stored since only outflow conditions are needed. 

In the same manner that shocks arise from the intersection of 
characteristic waves, they can also meet with other shocks to form new 
shocks. In addition, shocks introduced in the shock-forming zone of a 
given cascade element will propagate into the downstream elements of the 
cascade interacting with each other and creating new shocks. These shock 
interactions cannot be treated by the above PSF scheme, because it tracks 
only one wave at a time. For this reason, a further simplification is 
introduced which consists of restricting the shock interactions to the 
junctions of the kinematic cascade. The shocks emanating from the 
discretized inflow hydrograph are tracked across the shock forming zone 
using the explicit scheme, Eqns. 40. When all the shocks have been 
projected to the downstream boundary, their fronts are smoothed out using 
Eq. 41. A smooth outflow hydrograph, incorporating the overall effect of 
the shocks formed upstream is thus obtained. The interaction of shocks 
carried by converging outflow hydrographs is then simulated by simple 
superposition of these hydrographs. The resulting outflow hydrograph, 
which satisfies flow continuity, is in turn used as upstream boundary 
condition to the next element and the procedure is repeated. This method 
will be called an approximate shock fitting (ASF) scheme. This scheme 
computes the outflow hydrograph at any junction of the cascade reflecting 
the overall effects of the shocks formed in the upstream elements. 

2.4 SEDIMENT ROUTING 

Sediment yield from agricultural catchments is controlled by physical 
principles governing the detachment and transport of sediment particles. 
The source of sediment erosion, or sediment supply, is the detachment by 
raindrop impact and by runoff. If the sediment load from upstream areas is 
less than the potential transport capacity of the flow, the supply is 
depleted and erosion occurs. If the sediment load is greater, deposition 
occurs. These processes are all interrelated and must satisfy, locally, 
the conservation principle of sediment mass. Therefore, in addition to the 
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equations of continuity and momentum for water, additional equations are 
required for sediment routing. These are the sediment continuity equation, 
and relations describing sediment supply and transport. These equations 
are presented below, and they are equally applied to overland and channel 
flow. 


2.4.1 Sediment Continuity Equation 

The continuity equation for the size classes k = 1, 2, ..., n, forming 
the sediment load can be written as 


3G 


9x 


. 3C. A 

Sk + + (1-A) 


3Pz. 
_ V 

3t 


= g 


£k 


(42) 


where G . is the sediment transport rate by volume per unit time, X is the 

S K 

soil porosity, P is the wetted perimeter, z^ is the land surface elevation, 
g^k i s the lateral inflow, and 

C k = G sk /Q <M> 

is the fraction concentration by volume (Fig. 8a). The concentration and 
bed elevation for the total load are 


and 
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the third term in Eq. 42 will be denoted by 
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This term can be envisioned as a sediment supply term for exchange between 
the flow and the detached bed material during erosion or deposition. 
Assuming that the sediment moves essentially at the same average velocity 
of flow, V, Eq. 42 can be expressed as 
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Taking V as approximately constant over small space and time intervals, 
yields 


3A . 3A . 
sk , ,, sk 

3t V 3x 8 £k 8 dk 


(48) 


in which A , = C, A = G ./V is the volume of sediment fraction present in 

sk k sk 

the flow per unit length. Eq. 48 is used to track the aggradation or 
degradation of the bed as explained in section 2.4.4. 

The present model uses a power function to relate the wetted perimeter 
to the flow area. That is 

P = aA b 

where a and b are coefficients depending on the cross-sectional shape of 
the flow reach. In particular a = 1.0 and b = 0.0 for overland segments. 

2.4.2 S ediment Tra n sport Formu las 

These formulas are used to determine the sediment transport capacity 
for a specific set of flow and sediment characteristics. The formulas used 
in the present model were selected after a recent study by Alonso et al. 
(1981). They compared the predictions of nine transport formulas with 
flume and field data. The comparison was based on 40 field measurements, 
523 flume experiments, and 176 tests on concave slopes, with sediments 
ranging from coarse sands to very fine soil particles. As expected, no 
formula satisfactorily represented the entire spectrum of sediment and flow 
characteristics. Nevertheless, three of the tested formulas gave 

satisfactory estimates of transport capacity over different subsets of the 
data range. 

The total load formula of Yang (1973) best estimated streamflow 

carrying capacity in the range of fine to coarse sands. The total load 
formula of Laursen (1958) predicted reasonably well small streamflows 
carrying very fine sands and silts. It should be used with some 

reservations, however, for computing transport of lighter materials such as 
soil aggregates. The bed load formula of Yalin (1963) can be used to 

compute sediment transport capacities for overland flows. These 
conclusions are graphically summarized in Fig. 9. Each of these formulas 
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is presented below as the originator intended, but rewritten in terms of 
dimensionless parameters, (Alonso et al., 1981). Where the formulations 
required graphical solutions (i.e., determination of threshold conditions 
from Shield's curve), analytical equivalents (not shown here) have been 
worked out to facilitate their use in digital computation. 

(i) Total load formula of Yang (1973): Yang based his formula on the 
premise that total load is dominated by the rate of potential energy 
expenditure per unit weight of water. He used this concept and dimensional 
analysis to derive his formula, the coefficients of which were determined 
from a large set of laboratory data. This formula is: 

“’k = 6 k Z k (V/U * )(10<1> " 6/S] ’ (49) 

where 

«J» = 5.435 - 0.286 log (w k d k /v)-0.457 log (u*/w k ) + 

(1.799 - 0.409 log (w k d k /v) - 0.314 log (u*/w k )] 

log (VS Q /w k - V c S Q /w k ), (50) 

V c /w = 2.5/[log(u*d k /v) - 0.06] + 0.06, 0<u*d R /v<70, (51a) 

V /w = 2.05, u.d./v > 70. (51b) 

c w k 

(ii) Total Load Formula of Laursen (1958): based mostly on heuristic 
considerations, Laursen developed a formula relating the load concentration 
to the relative roughness and excess tractive force. This relationship is 
corrected by an empirical function of (uj,./w) which accounts for the 
effectiveness of turbulence in suspending the bed material. The 
contributions of each size fraction are added up to yield the total load. 
The formula is: 
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Fig. 9 Range of applicability of sediment transport formulas 











where f(u^/w^) is the empirical function obtained by Laurstn from t 1 um< 
experiments. 

(iii) Bed Load Function of Yalin (1963): this formula is based on a 
theoretical analysis of saltating particles, in which it is proposed that 
the bed load rate is related to the range of the particles' saltation 
rather than to the number of particles participating in the motion fhr 
resulting formula is of the excess-shear type, and was calibrated on a 
limited set of laboratory data. This formula has the form 

* k = 0.635 0 1' 2 II - (0 ck /e k )l |(i/e ck ) 

- (1/ao 0 cR ) ln( 1 + ao] , (5.3) 


in which 


a = 2.45 S" 2/S 0^. , (54) 

ck 

a = (e 5 O /0 ck ) _ 1 = fu */e c k (s ‘i) gd k ] ' 

In Eqns. 49 through 55 <t> k is the dimensionless volume transport rate, 
0^ and are the mobility number and relative roughness based on the d^ 
fraction size, w^ is the settling velocity of sediment, and the subscript c 
denotes threshold conditions. These parameters are defined as 


\ = (Gsk/yPSi/Us-Dgdj^, 

( r .b) 

6 k = f C s -i)g d k ], 

(57) 

Z k = y/d k> 

(58.) 


where S is the specific gravity of the sediment fraction, u. v is the bed 
shear velocity, and y and V are the specific weight and kinematic viscosity 
of water, respectively. 
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2.4.3 Equations of Sediment Supply 

Soil eroded from land areas with distributed and concentrated flows is 
the source of most of the sediment transported in a catchment system. Soil 
erosion is a complex process of detaching soil particles and transportng 
them downslope through the action of raindrop impact and runoff. Erosion 
begins when raindrops strike the land surface and detach soil particles by 
splash. Whenever the soil surface is not protected by vegetation, or any 
other cover, raindrops can detach very large quantities of soil. Most of 
this eroded soil is moved downstream by surface runoff. When runoff 
reaches sufficient intensity, the rate of soil erosion becomes also 
dependent on the runoff characteristics and on the susceptibility of the 
soil to the forces of the flowing water. The following sections describe 
the approaches used in the present model to simulate these major sources of 
eroded sediment. 


2.4.3.1 Soil Detachment by Raindrop Impact: Past research on the process 
of soil detachment by raindrop impact suggested that the rate of detachment 
is proportional to the square of the net rainfall intensity (Meyer and 
Wischmeier, 1969). However, bare soils differ in their susceptibility to 
detachment by raindrop impact due to a variety of properties such as 
primary particle size distribution, soil structure, organic matter content, 
etc. Therefore, an erodibility parameter, a^, is introduced to describe 
the rate of erosion as 

D = a I 2 . (59) 
r r 

This relationship has been recently cor.i.rmed by Meyer (1980) who reported 
extensive field studies showing that Eq. 59 is a valid approximation for a 
wide range of soils and cropping conditions. Eq. 59 is the basic equation 
for rainfall detachment. Parameters must be added to account for other 
factors that influence the rate of detachment. Meyer (1980) has shown that 
the I 2 law is not affected by the stage of canopy development, but the a f 
coefficient is dependent on the erodibility of the soil and decreases as 
the vegetative cover changes from first growth to full canopy. To account 
for this effect a cover factor is added based on the cropping-management 
factor C used in the USLE (Wischmeier and Smith, 1978). The factor C is 
defined as the ratio of soil loss from land cropped under specified 
conditions to the corresponding loss from tilted, continuous fallow. 











Wischmeier (1972) presented a method and the necessary graphs for 

determining C for situations were data are not readily available. In this 

method C is treated as the product of three subfactors depending on (i) 

canopy cover, (ii) mulch and ground cover, and (iii) residual effects of 

the land use, Cjjj, respectively. Approximating the first two subfactors 

as linear functions of the densities D and D introduced in section 2.1, 

c g * 

the factor C can be expressed as 

C = 0.2C in ((l-D c )(l-D g ) + H c D c (l-D g )] 

where is a function of the average fall height of drops from the canopy 

cover. Introducing C in Eq. 59, neglecting the term depending on H for 

simplicity, and incorporating the coefficient 0.20^ into a^, yields 

D = a I* (1-D )(1-D ). (60) 

r r n c g 

Ponded water deeper than a critical depth cushions the impact of raindrops 
and also diminishes the erosion caused by the dissipation of impact energy. 
Mutchler and Young (1975) suggested that a water depth of more than three 
times the median drop size essentially eliminated detachment by raindrop 
impact. Consequently, detachment can take place only if the raindrops can 
penetrate through the water depth, h, and the thickness of detached soil, 
e, accumulated from previous events. Laws and Parsons (1943) found that 
the median drop size, expressed in millimeters, is related to rainfall 
intensity as 


°50 = 2.23 I 


0.182 


(61) 


Using this formula, the effect of water pondage is described by modifying 
Eq. 60 as 


D r = a r *n U'^) (l"D g ) I l-(h + e)/3D 5Q ], if 

(h + e) g 3D 50 , (62) 

and D r = 0 otherwise. A similar expression has been proposed by Li (1979). 
The amount of soil detached by raindrop impact during a time step At, and 
to be added to the current storage of detached soil is, therefore, 



Ae rk = a r IJ(l-D c )(l-D )(l-(h+e)/3D 50 Jp k At, (h+e) S 3D 5Q , (63a) 

Ae rk = 0, (h+e) > 3D 5q , (63b) 

where p k is the percentage of sediment material in the kth-size fraction. 

2.4.3.2 Soil Detachment by Flow: Erosion by overland flows usually 

occurs in many single rills. However, for modeling applications rill 

erosion is assumed to be uniformly distributed overland, and erosion by 

concentrated flows is restricted to channels. Erosion by flow potentially 

occurs if the sediment load carried by the incoming flow is less than the 

transport capacity of the flow. If the sediment load is greater, 

deposition occurs. Erosion by flow is assumed to occur at capacity rate if 

no sediment is present in the flow. But if the transport capacity of the 

flow is partially filled a corresponding depletion of the detached soil 

available for transport is computed. If the transport capacity is less 

than the available detached soil, transport is the limiting factor and no 

erosion by flow takes place. If the available detached soil is less than 

the transport capacity, additional soil is detached by the flow to 

compensate for the insufficient supply. These concepts are implemented as 

follows. The amount of sediment that the flow can potentially carry 

downstream during the time increment At is transformed into an equivalent 
p 

thickness Ae , as explained in the following section. If this thickness is 

less than the layer of available detached soil, e (Fig. 8a), no detachment 

p 

by flow occurs. If Ae > e, the available detached soil is not enough for 

transport and the flow can potentially detach the additional amount e - 
p 

Ae . However, the resistance of the soil to the erosive forces of the flow 
depends on the soil properties and structure, as well as on the condition 
of the land surface (Olson and Wischmeier, 1963; Kemper, 1966; Wischmeier 
and Mannering, 1969; Grissinger, 1980). Consequently, the potential 
thickness of detachment by flow is modified by a flow detachment 
coefficient, a*, yielding the following amount of soil detached by flow: 










(64) 


**fk = * a f (AeP ' e) Pk* 

where ranges from zero to 1.0 depending on the soil erodibility. The 
detachment coefficients in Eqs. 63 and 64 are optimization parameters that 
are calibrated by fitting the sediment discharge rates to observed data. 


2.4.4 Numerical Procedure for Sediment Routing 

The equations of sediment supply and sediment mass conservation are 
coupled to the water routing scheme using the procedure described below to 
track the evolution of the sediment load, and to compute the aggradation 
and degradation of the land surface and stream channels. 

Eq. 48 is a linear hyperbolic equation that may be solved by the 
method of characteristics. The steps given below follow those presented in 
section 2-3.2. The total differential of A gk (x,t), 


3A 


sk 


8x 


8A . 

dx + ~0t dt ~ dA sk , 


(65) 


and Eq. 48 form a system of two equations in the two unknown partial 
derivatives of A gk> that is 

"l V 
dt dx 

Since these derivatives do not exist along the characteristic paths, the 
determinant of the coefficient matrix vanishes at points on the 
characteristics. This condition yields the characteristic equation 


8 \k /8t 
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( 66 ) 


Tt = V * Q'A 


(67) 


in agreement with the assumption that the sediment moves with the same 
velocity of the flow. One of the invariants of this solution also gives 


dA 


sk 


dt 


= g £k + g dk 


( 68 ) 
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or, by integrating Eq. 68 


sk 


(x ’t) = A sk (x o’ t o ) + { (g 
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«dk>«- 


(69) 


where A sk (x o >t o ) is an initial value of Applying this equation to 
two points E(x^_j,tj_j) and F(x^,t^) on a characteristic path (Fig. 9) and 
expressing the result in discrete form, gives 


^ A sk\,j ^ A sk^i-l,j-l + 


“Wi.j * 




(70) 


in which g^ and g^ are assumed uniformly distribured in the interval x^ 

p 

= x 5 x^. The potential transport capacity of the flow, C^, is obtained 
from the formulas in section 2.4.2 using the average of the flow 
characteristics computed at E and F, and the sediment properties at point 

p - 

E. Substituting (C^ A) for the right-hand-side of Eq. 70 and solving for 

^dk^.j yields 


(g dk^i,j " At. fC k A * ^ A sk^i-1, j-1 ^ ' ^Ak^i, 


(71) 


p 

where the overbar denotes the average flow area, and g^ is the potential 

change in detached soil storage caused by either erosion or deposition, 
p 

Deposition : If < 0, the potential carrying capacity of the flow is 

less than the sediment present in the flow. Consequently, the sediment 
load transport rate is 


Q (V, j = Q < . (72) 

P 

and the excess load, - g^i is deposited on the bed. However, whether this 
amount will reach the bed during the time step At depends on this being not 
less than the average time for the sediment particles to settle to the 
channel bed. Data by Jobson and Sayre (1970) and by Lean (1971) indicate 
that this settling time may be computed using the particle fall velocities 
in the quiescent fluid. Therefore, the actual deposition of a size 
fraction during the interval At is calculated from 
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or 


«Vi,j * ■«**'. 


‘Vi.j * 'KWij • lf * < '• 


where P = w^ At/h. The new bed elevation is 
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n <V. 


Z. . = Z. . . + - l . 

i.J i»J"l p k=1 1-* 


ill 


(73a) 

(73b) 


(74) 


p 

Erosion : If > 0 the potential carrying capacity exceeds the amount of 

material in transport and, therefore, the flow will entrain additional 
material. In this case, one of two possible events may occur depending on 

p 

whether g^ is equal to or greater than the thickness, (e^^ _j, of the 

available detached soil. Let 


(e k } i,j = < ’ e k ) i,j-1 + ^rk^i.j-l ' 


(75) 


denote the depth of detached soil in storage at the end of the time step 

At., in which 
J 


K>i,j = <4\,j v f 


(76) 


is the depth of potential erosion by flow. 

(i) If (e. ). . £ 0 the available detached soil is enough to supply 

k i , j 

the sediment entrained by the flow. In this case no erosion (of undetached 
soil) occurs, the transport rate at the end of the time step At^ is equal 
to the carrying capacity (Eq. 72), and (Z, ). . = (Z.). . . 

(ii) If (e, ). . < 0 the availability of detached soil is less than 

k i , j 

the potential entrainment, and additional soil is detached by the flow. 
The depth of erosion by flow is 


(Ae,,). . 
fk'i.j 


= -V'k’i.j 


(77) 
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from which 


'Wi.j = + 


(78) 


From Eqns. 70 and 78, the concentration of the updated load is 


(A ,). . At. 

(CJ; 4 = + -T 1 (8.J,- 
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(79) 


Finally, the new elevation of the eroded bed is 


n ^ e fk^i i 
Z. . = Z. . , + I - 

k=l l-\ 


(80) 


Substituting Eq. 72 or 79 into Eq. 44 gives the concentration of the total 
sediment load. The total concentrations existing on all the 
characteristics paths reaching a certain location define the outflow 
sedimentgraph at that location. 

In some instances, the time step size selected for water routing may 
yield a cluster of small space increments (Fig. 8b). Repeating the 
sediment calculations for each of these increments may not be necessary 
when the sediment routing parameters are only required at points spaced 
farther apart. In these cases it is more efficient to combine a number of 
space and time steps used for water routing, into larger steps for sediment 
routing. To this end, the program contains an option to compute the 
sediment routing parameters only at points satisfying the conditions 
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where GDX is a user supplied parameter (Fig. 8b). 
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(81b) 
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APPLICATIONS 


To illustrate the relative accuracy of the present model, it is 
applied to several examples reported in the literature. These are the 
hypothetical three-plane cascade discussed by Kibler and Woolhiser (1970); 
the experiments performed by Iwagaki (1955) involving unsteady, open- 
channel flow with lateral inflow; the runoff studies carried out by Schaake 
(1970) in a small urban catchment; and two agricultural catchments reported 
by Burford and Clark (1973). The PSF and ASF schemes are compared on the 
three plane cascade example. The results from the ASF shock-fitting scheme 
are also compared with those obtained with the implicit finite-difference 
(FD) scheme presented by Li et al. (1975a). 

The computational parameters used in each example are given in Table 
1. The kinematic wave parameter n, Eq. 26, is kept fixed at 1.50 (Kibler 
and Woolhiser, 1970; Singh, 1975). With n fixed, the calculations only 
required characterization of the parameter a. The values of a, estimated 
by Kibler and Woolhiser (1970), are used in the three-plane cascade 
computations. In all the other examples, this parameter is computed, using 
the relationship proposed by Singh (1975) 

a = Cj + C 2 (S Q )* , (82) 

where C^ and C^ are constants to be optimized in each case. This 
relationship was chosen only for its simplicity. However, any other 
expression could have been used, since the computational efficiency of the 
schemes is not contingent on the manner or is estimated. 

3.1 HYPOTHETICAL KINEMATIC CASCADE 

Kibler and Woolhiser (1970) used a three-plane kinematic cascade to 
illustrate the formation of kinematic shock waves. The planes are 400 ft. 
long with slopes 0.04, 0.01, and 0.0025, respectively. The lateral inflow 
is a rainfall pulse with an intensity of 0.75 inch/hr and a duration of 30 
min. The necessary condition (Eq. 37) is satisfied at the unions of the 
three planes. Thus, shock-forming regions exist in the rising limbs of the 
upstream inflows to the second and third planes. 

The characteristics and shock paths crossing the x-t domain of the 
cascade are presented in Fig. 10. The shock paths were computed using the 
PSF scheme. Because there is no upstream inflow to the first plane, the 
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outflow from this plane is a smooth hydrograph generated by the 
characteristic waves emanating from dry ground conditions. This hydrograph 
is used as inflow to the second plane. The rising limb satisfies the shock 
forming conditions (Eq. 36), and so, small shocks are introduced by 
discretizing this portion of the hydrograph. These shocks are projected to 
the downstream boundary along with the characteristics emanating from the 
initial dry ground condition. There they define an outflow hydrograph 
containing a smooth rising limb formed by the characteristic waves, 
followed by a piecewise continuous part formed by the shock waves. To 
reduce the number of shocks to be traced in this example, shocks arriving 
within the same time step are assumed to merge and proceed as a single 
shock to the next plane. The same procedure is applied in routing over the 
third plane. The hydrograph at the cascade outlet contains all the shocks 
crossing the last two planes. The shocks slow down as they enter the third 
plane, reflecting the decrease in the value of a. The complete outflow 
hydrograph is compared in Fig. 11 to the analytical solution of Kibler and 
Woolhiser (1970). Their solution displays the two shocks emanating from 
the unions of the planes. The two hydrographs are quite close, which shows 
that the approximation introduced in formulating the PSF solution does not 
detract from its accuracy. 

The same hydrograph calculation was carried out using the ASF and FD 
schemes. The ASF solution is plotted in Fig. 11. It exhibits a smooth 
rising limb showing the overall effect of the shock waves, and is in good 
agreement with the two other solutions. Moreover, it is twice as fast as 
the PSF scheme (Table 1). It is thus evident that the ASF scheme, in 
addition to being easier to work with, gives results as accurate as the 
analytical solution. For these reasons the PSF technique was abandoned in 
favor of the more efficient ASF solution. 

The hydrograph computed, using the FD scheme with the same values of 
a, is compared in Fig. 12 with the ASF solution. Although the hydrographs 
agree in their predicted yields, the FD solution exhibits a pronounced 
delay and reduction of the peak. In an attempt to improve this solution, 
the FD calculation was repeated using a new optimized value of the 
kinematic parameter. The result shows improvement in the peak estimate but 
an overall deterioration of the hydrograph shape (Fig. 12). This inability 
of the FD scheme to reproduce the analytical solution is further discussed 
in the context of the next example. 
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Fig. 12 Comparison of finite-difference results with approximate 
shock-fitting solution 










3.2 


UNSTEADY CHANNEL FLOW 


Iwagaki (1955) reported a set of laboratory experiments of unsteady 
flow with lateral inflow in a smooth open channel 7.3 ft long. The channel 
was divided into three sections of equal lengths and different slopes. 
From the upstream end, the slopes were 0.020, 0.015, and 0.010. The 
lateral inflow was adjusted so that each section would receive a constant 
rate of inflow, with the upper, middle, and lower receiving 0.0425, 0.0251, 
and 0.0315 inch/sec. The durations of lateral inflow used in the 
experiments were 10, 20, and 30 sec. Under these flow conditions, Eq. 37 
predicts a shock forming zone at the union of the upper and middle 
sections. 

Figs. 13, 14, and 15 show that the hydrographs computed with the ASF 
scheme are in good agreements with the partial-equilibrium hydrographs, 
measured by Iwagaki. The kinematic parameter a was adjusted by fitting the 
data plotted in Fig. 13. Because the slopes, geometries, and hydraulic 
roughness were the same in all three experiments, the same value of a was 
used in the calculation of the other two hydrographs. These solutions 
correctly simulated the shock formation and its arrival time at the outlet. 
This is clearly depicted in Fig. 15, where the hydrograph peaks immediately 
after the lateral inflow terminates and falls continuously, until a sudden 
rise occurs when the shock arrives. Fig. 14 shows the shock arriving 
shortly after the first peak, whereas in Fig. 13 the shock arrives well 
ahead of the peak. 

Two different hydrographs computed with the FD scheme are given in the 
same figures. The hydrographs computed with the values of a used ii. the 
ASF solution do not agree very well with the measurements and, in addition, 
they do not reproduce the aforementioned shock effect. A second solution 
was obtained by recalibrating the kinematic parameter but the new 
hydrographs do not exhibit any better agreement. In the present example, 
the shock formation is an intrinsic part of the solution to the kinematic 
wave approximation. However, because of its smoothing effect, the FD 
scheme is unable to reproduce this important aspect of the analytical 
solution. In addition, this scheme required between 50 and 100 percent 
more computing time that the ASF solution (Table 1). 
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13 Simulation of Iwagaki's partial-equilibrium hydrograph for lateral 
inflow duration of 30 sec 











Fig. 14 Simulation of Iwagaki's partial-equi1ibrium hydrograph for lateral 
inflow duration of 20 sec 













ig. 15 Simulation of Iwagaki's partial-equilibrium hydrograph for lateral inflow duration of 10 sec 












3.3 


URBAN CATCHMENT 


This example is used to illustrate that the shock-fitting technique 
presented in this report is applicable to a combination of interconnected 
overland and channel segments. The case considered herein is the small 
urban catchment, reported by Schaake (1970). In this publication Schaake 
describes an event labeled 3 SPL1 on a 0.39 acre impervious parking lot 
labeled SPL1. He represented the catchment by a number of interconnected 
overland segments and V shaped channels (Fig. 16). The overland segments 
vary from 20 to 36 ft in length, and their slopes range from 0.0167 to 
0.019. The V-shaped channels have lengths varying from 50 to 165 ft, and 
slopes ranging from 0.0148 to 0.0213. The channel-side slopes were all 
1:11.3. The same geometric representation is used in the simulations 
reported herein. 

The runoff hydrographs measured at the outlet of the catchment and the 
hydrographs computed with the ASF and FD schemes are shown in Fig. 17. The 
value of a was obtained by fitting the measured data. The result obtained 
with the ASF scheme overpredicts the first peak but agrees very well with 
the rest of the data. The FD calculation obtained with the same a does not 
predict the second peak well. In an attempt to improve the second peak, a 
second calculation using the FD scheme was made by recalibrating a. This 
improved the second peak but did not help the rest of the hydrograph (Fig. 

17) . In general, the hydrographs predicted by the FD scheme do not exhibit 
the significant deviations shown in the previous calculations. This is 
because, in this example the lateral inflow rate changes rather smoothly, 
thus reducing the magnitudes and effects of shocks. Also, the slopes in 
the previous two examples exaggerate the development of shocks, whereas in 
this example the slopes are nearly the same, thus leading to much smaller 
shocks. Nevertheless, these solutions take twice as much computing time as 
the ASF calculation (Table 1). 

3.4 AGRICULTURAL CATCHMENTS 

The data used in these tests were obtained from the USDA experimental 
catchments W-5 southwest of Holly Springs, MS, and R-5 near Chickasha, OK 
(Burford and Clark, 1973). Catchment W-5 drains a 1.76 mile 2 area (Fig. 

18) , with a good mixturp of cultivated land, timber, pasture, and idle 
land. Catchment R-5 has an area of 23.7 acres, and is range land with an 
excellent native grass cover (Fig. 19). These catchments were chosen 
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because of diversity and availability of most of the required data. 
Sediment records were available for catchment R-5, but the sediment yield 
was so small that no comparisons were made between measured and computed 
sediment discharges. 

The model was calibrated using the event of February 21, 1971, on 
catchment W-5, and the event of May 6, 1969, on catchment R-5. 

Infiltration parameters for catchment W-5 were estimated from information 
reported by Smith and Parlange (1978) for Colby swelling type soils, 
because very little infiltration data were available. The infiltration 
parameters for catchment R-5 were obtained from an average infiltration 
curve obtained from a large number of infiltrometer runs conducted in the 
fall of 1977. The calibrations were carried out by adjusting first the 
flow resistance parameters (Eq. 82) to match the hydrographs, and then the 
sediment-model parameters (Eqns. 63 and 64) were adjusted to fit the 

sedimentgraphs. These same parameters were used in simulating all the 
remaining events. The results of these tests were reported in an earlier 
paper by Alonso et al. (1978). 

Examples of the comparison between the simulated and the measured 
hydrographs and sedimentgraphs are shown in Figs. 20, 21, and 22. A total 
of nine events on catchment W-5 and two on catchment R-5 were simulated. 
The agreement between the shapes of measured and simulated events is 
satisfactory. Comparisons between measured and computed water yield, peak 
runoff rates, sediment yield, and peak sediment discharge are given in 
Figs. 23 and 24. These plots show that the model estimates both yields and 
peaks within a range of about ±40 percent of measured values. The limited 
number of events used in this study precluded an estimation of the 
confidence level of the range of variability. The results presented in 
Figs. 20 through 24 indicate that simulations of different size events 
using only one set of parameters for each catchment, were satisfactory. 
This suggests that the model could be used to predict the response of a 
catchment to different management practices, if the model parameters 
associated with each practice could be accurately estimated. Also, the 
above results indicate that the model could be readily transferred to 
ungaged catchments, if the model parameters were properly regionalized. 
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Fig. 21 Sedimentgraph from catchment W-S for the February 21, 197] event 
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CONCLUSIONS AND RECOMMENDATIONS 


4 

4.1 CONCLUSIONS 

1. A numerical model for routing water and sediment on small 
catchments has been developed. The model accepts any single rainfall 
hyetograph and produces runoff and sediment hydrographs for the modeled 
catchment. 

2. The model is developed on a general basis so that it may be 
applied to any agricultural catchment by changing only the input data. The 
approximate range of basins over which the model is applicable is from a 
few acres to about 5 square miles. 

3. The model is based on the physical processes governing the 
mechanics of water and sediment movement and requires the calibration of 
four parameters. 

4. The model can be used to simulate the effect of different land 
uses on the water and sediment yields from the modeled catchment. 

5. The model predicts the surface component only. It does not 
presently predict subsurface and groundwater movement. 

6. The applicability of the model is restricted to streams where the 
channel geometry does not change significantly during a storm event, and in 
which the kinematic-wave approximation for flow routing is valid. 

7. The present model simulates single storm events; the user must 
estimate the initial conditions for the storm. The model can still be 
applied to a sequence of rainfall events if the user can make satisfactory 
estimates of the initial soil moisture conditions. In this case the model 
can be used to predict a sequence of surface runoff and sediment transport 
events. 

8. The model has been validated with several data sets including data 
from the natural catchments W-5 in northern Mississippi, and R-5 near 
Chickasha, Oklahoma. The shape of water and sediment hydrograph and total 
water and sediment yields of a number of events were satisfactorily 
simulated. 

4.2 RECOMMENDATIONS 

1. It is recommended that the model be put to work on real systems. 
The evaluation and continuous updating of the model are essential to its 
credibility and effectiveness. 
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2. It is recommended that the model be further developed, or 
modified, to permit continuous simulation over long time spans (20 to 50 
years). This is essential in evaluating the long term response of a 
catchment, or stream network, which is dependent not only on the history of 
management practices, but also on the sequence of storm events. 

3. It is recommended that the model be further developed to track the 
channel geometry of streams that become unstable due to bank erosion and 
deposition. 

4. It is recommended that data gathering efforts be continued to 
provide an adequate base for further model development and validation. 
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ADDENDUM 1. DESCRIPTION OF DATA INPUT AND COMPUTER PROGRAM 


The structure of the program SEDLAB, a software system developed for 
the simulation of water and sediment movements in agricultural catchments 
is described. The numerical schemes on which the program is based have 
been described in detail in Part 2. A description of the data input and 
important variables used in the program is given in the following sections 
of this addendum. The last section presents a complete list of the 
program. 

The system SEDLAB consists of a main program and several subroutines. 
The sequence of operations performed by the main program is shown in Fig. 
1.1. The main program inputs the required data to the system, calls the 
subroutines according to the computation scheme, and prints out the 
calculated results. Subroutine INTCPN computes the evaporation and 
interception losses, and determines the net rainfall rate. Subroutine 
INFLTN computes the ponding time, and the infiltration losses in a segment 
for one time step. Subroutine WROUT routes flow through a segment for one 
time step and it calls in turn subroutine SROUT. This subroutine performs 
sediment routing through a segment for one time step. Subroutines LAURSN, 
SETVEL, SHIELD, YALIN, and YANG are called by SROUT to compute potential 
carrying capacities and sediment transport parameters. The program has 
been designed to process the entire network of segments for one time step 
according to the computational sequence described in section 1.2 of this 
addendum. Once the entire network has been processed the simulation neves 
to the next time step. This sequence is repeated until the entire 
simulation period is computed. 

The computer code requires 49,710 words on a Mod Comp Classic computer 
system. This is a 16-bit machine that uses two words for each single 
precision variable. The execution time for the event of February 21, 1971, 
on catchm> v-5 is 90 seconds. 

A! 1 DATA INPUT 

'uta inquired to run program SEDLAB includes: 

«• !(• i i c.t, length, bed slope, bed elevation, and wetted 

• . i. , • i- * i »■ 1 at i ons . 

. • . if. ( . and ground cover density, interception 
•• ! ground (overs, and ratio of evaporating 
, - I f : g 11 und ruver. 























Soil data : soil saturated hydraulic conductivity, soil sorptivity, 
specific gravity and si 2 e distribution of bed material. 

Flow and sediment routing parameters : constants describing flow 

resistance, parameters describing soil detachment by raindrop impact and 
surface runoff, maximum penetration depth of raindrop impact, and 
computational sequence. 

Storm characteristics : rainfall intensity, mean evaporation rate, and 
initial interception storage content. 

A detailed description of the data input is given below: 


Card FORTRAN 

No. Variable 

Descript ion 

Units 

Format 

1 TITLE 

Alphanumerical identification of 

simulation run. 

- 

20A4 

2 AREA 

Drainage area of catchment 

acres 

F10.4 

NOV 

Total number of overland segments 

- 

14 

NCH 

Total number of channel segments 

with negligible infiltration 

"" 

14 

NCI 

Total number of channel segments 

with significant infiltration 

— 

14 

NSTHM 

Total number of storm events 

simulated in the run 


14 

3 (This card must he repeated lor each overland 

segment) 


SEG 

Number identliving overland segment 

Segments are numbered sequentially 

from 1 to NOV. 


14 

OVA 

Area of overland segment 

ft* 

F12.3 

SEEN 

Slope length of overland segment. 

This length is computed by dividing 

OVA by the length of the receiving 

channel reach. 

ft 

F10.4 

SLOPE 

Slope of overland segment. 

ft/ft 

F10.4 

CPER 

Coefficient in the wetted perimeter 

versus flow area relation (the 

default value is one) 


F10.4 


I .83 









EPER 


F10.4 


Exponent in the wetted perimeter 
versus flow area relation (the 

default value is zero) _ 

(This card must be repeated for each channel segment) 

SEG Number identifying channel segments. - 14 

These segments are numbered sequentially 
from (NOV+1) to (NOV+NCH) if there are 
no channel reaches with significant 
infiltration. Otherwise, the channel 
segments are numbered from (NOV+1) to 
(NOV+NCH+NCI). 


SLEN 

SLOPE 

CPER 

EPER 

Length of channel segment 

Bed slope of channel segment 

Coefficient in the wetter perimeter 

versus flow area relation. 

Exponent in the wetted perimeter 

versus flow area relation. 

ft 

ft/ft 

F10.4 

F10.4 

F10.4 

F10.4 

VOG 

Interception storage capacity for a 

typical ground cover. 

inches 

F10.4 

SRG 

Ratio of evaporating surface to the 

horizontal projected area of typical 

ground cover 

it 1 lit 2 

F10.4 

VOR 

Ratio of the interception storage 

capacity of a typical canopy cover to 

that of a typical ground cover. 


F10.4 

HLR 

Average height of ground cover in 

channels. 

ft 

F10.4 

NSOIL 

Number of representative sediment 

size fractions used in the simulation. 

- 

14 

SPGR 

Specific gravity of sediment 

- 

F10.4 

AIM 

Coefficient of soil detachment by 

- 

F10.8 


raindrop impact (a^). User supplied 
optimization parameter. 


ADF Coefficient of soil erosion by surface 

flow (a^). User supplied optimization 
parameter. 


F10.8 












GMAX 

Maximum penetration depth of 

raindrops (Eq. 61). 

ft 

F10.4 


GDX 

Space increment adopted for sediment 

routing. User supplied parameter. 

ft 

F10.4 

7 

D50 

Median size of sediment fractions. 

mm 

10F7.3 



The program can accomodate up to 





five fractions. 



8 

PC 

Percentage of sediment fractions. 

* 

10F7.3 

9 

(This 

card must be repeated for each segment in 

the sequences used 


for ' 

cards 3 and 4) 




CANO 

Canopy cover density for the segment 

ft 2 /ft 2 

F10.4 


GCOV 

Ground cover density for the segment 

ft 2 /ft 2 

F10.4 


HYCND 

Saturated hydraulic conductivity for 

the segment 

in/hr. 

F10.4 

10 

(This 

card must be repeated for eaih segment) 




ISEG 

Index identifying the position of the 

segment in the computational sequence 

(see following section) 


14 


I ARY 

Array containing the storage locations 

- 

514 



of the inflows and outflows computed lor 




the segment (see following section) 



11 

TEMP 

Water temperature 

Celsius 

F10.4 


GAMA 

Specifu weight uf w.itci 

lb/ 11 3 

no.4 


SNU 

Kinemalu visiosity <>t water 

(t 2 /sec 

F10.8 


EXP 

Exponent in fq. 2<> 

- 

no.4 


Cl 

First eoel lie lent desinhing 

- 

F10.4 



ki Hemal n -wave f rut ion parameter 





(F.q. 82 1. User supplied optimization 





pa rameter 




C2 

Second coefficient in Eq. 82. 

User supplied optimization parameter 

- 

F10.4 

12 

NEED(I) Vector representing the following 

- 

614 


output options: 
NEED (1) input data 
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NEED (2) Bed elevation changes computed for 
each segment at the end of the 
simulation period 

NEED (3) water budget and sediment yield 
NEED (4) Computed infiltration rates. Measured 

and computed hydrographs and sedimentgraphs 
NEED (5) Plots of hyetograph and hydrograph 
NEED (6) Plots of hyetograph and sedimentgraph 


When NEED(I) > 0, the program prints selected output; 
when NEED(I) = 0 the program does not print selected 
output. _ 


13 (Cards 

13 through 17 must be repeated for each 

storm event) 


STORM 

Alphanumerical identification of the 

storm event 


5A4 

DTM 

Size of time step 

min 

F10.4 

ITMAX 

Number of time steps in rainfall 

duration 


14 

ITCOM 

Number of time steps in simulation 

period 


14 

EVP 

Mean evaporation rate 

in/hr 

F10.4 

VIN 

Initial interception storage, defined 

in/in 

F10.4 


as the ratio of the initial storage 
to the intercept ion st orage c apacity 


14 

(This card must be repeated every eight segments until all 

segments are included) 

SORPTY Sorptivity parameter for each of the in 2 /hr 

eight consecutive segments (S z /2, 

Fig. 3) 

the 

8F7.7 

15 

(This card must be repeated every twelve time 

rainfall intensities are read in) 

steps until 

all the 


DR Rainfall intensity for each of the 

twelve consecutive time steps 

in/hr 

12F6.3 

16 

(This card must be repeated every twelve time 

steps are included) 

steps until 

all ITCOM 


QMES Water discharge measured at the 

catchment outlet for each of the 

twelve consecutive time steps 

ft 3 /sec 

12F6.2 


1.86 












17 (This card must be repeated every twelve time steps until all ITCOM 
steps are included) 

GMES Total sediment discharge measured at Ibs/sec 12F6.2 

the catchment outlet for each of the 

_ twelve consecutive time steps _ 

A1.2 PREPARATION OF INPUT ARRAYS ISEG AND IARY 

The first step in preparing these data is to divide the catchment into 
interconnected overland and channel segments, each homogeneous within 
itself, and to assign an identification number to each segment. The 
principal direction of flow is determined for each overland segment from 
the contour line map of the catchment. The flow path through the cascade 
of segments is established by following the logics of gravity and flow 
continuity. The order in which the segment numbers appear in the flow path 
defines the computational sequence. This procedure is illustrated in Fig. 
1.2 which presents the segmentation used in simulating the catchment W-S 
described in section 3.4. The arrows shown in the figure denote the flow 
direction in each overland segment. It should be noticed that the 
segmentation is restricted to no more than two overland flow segments as 
input to a channel segment, and at most two inflow channel segments to a 
downstream channel segment. 

After the segmentation of the catchment has been completed, two arrays 
must be 6et up by the user. The first, [SEG, tells the program the 
sequence for computing flow and sediment discharge for each segment. The 

other, IARY, tells the program for any segment where to find previously 

computed inflows and where to store the computed outflows. These inflows 
and outflows are stored in different columns of the matrices Q and GS 
described in the next section. 

A table, AUX (Fig. 1.3), is used as an aid to set up the two arrays. 
AUX and the matrices Q and GS have the same number of columns. The table 
is constructed as follows. Starting with one of the most upstream channel 
segments, its inflow segments are selected. Then their numbers are placed 
in the first available columns of AUX and in separate rows. This will 

usually mean an overland flow segment number in row 1, column 1, and 

another in row 2, column 2. The channel segment number is then placed in 
the next available row and the first available column if no further inflows 
to the segment need to be computed. An available column is one that does 
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1.2 Geometric segmentation and flow path for catchment W-5 















not have the number of a segment waiting to be combined in some further 
downstream segment. This procedure is continued through the flow path, 
inserting segment numbers in the table until the outlet of the catchment is 
reached. If a junction with another channel segment is reached, then what 
has been computed must be retained while the area upstream of the junction 
is computed. 

To illustrate consider Fig. 1.2. Channel segment 29 is chosen as the 
first most upstream channel segment. Other segments which could have been 
chosen as well to start are 30, 33, 34, 35, or 20. The overland flow 
segments for segment 29 are 1 and 2. These are placed respectively in row 
1, column 1, and row 2, column 2. When the two overland flow segments are 
combined to give the channel outflow in 29, they are no longer needed. 
This allows the outflow from 29 to be t tered in column 1 of the next row. 
Next, the outflow from channel segment 30 must be computed. This requires 
computing the overland flow segments 3 and 4. They are entered in columns 
2 and 3 because 29 must be saved and is in the column 1. When 30 is 
computed it must be saved also until the inputs 5 and 6 to segment 31 have 
been computed. When 29, 30, 5, and 6 have all been computed, then 31 can 
be computed and placed in the first column of the next row. None of the 
segments previously computed need to be retained at this point. When 
segment 32 has been computed it must be retained until segment 39 has been 
computed which requires moving to upstream segments above 39. The channel 
segments 12, 33, and 38, and their corresponding lateral inflows, are the 
upstream starting points above 39. 

In the given example, the segment numbers entered in AUX occupy five 
columns. However, a finer segmentation or a denser drainage network would 
have resulted in a larger number of columns. The number of rows of AUX is 
always equal to the total number of segments because each segment is 
computed only once in each time step. 

When AUX is completed, ISF.G and IARY are constructed from it and then 
AUX is discarded. To construct the column matrix ISEG, the number of the 
segment in each row of AUX is placed in the corresponding row of ISEG (Fig. 
1.3). The order of segment numbers in ISEG gives the computational 
sequence used by the program to route water and sediment through the flow 
path. 
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Fig. 1.3 Assemblage of input arrays ISEG and IARY for catchment W-5 
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To construct IARY, each segment in the catchment is considered. The 
information for segment number n is placed in the nth-row of IARY. Columns 
1 and 2 of IARY contain information about any overland flow segments which 
are input to that segment. These columns will contain a zero if there are 
no overland flow segments, and will contain the column locations for the 
overland flow segments in AUX if there are (Fig. 1.3). For example, 
segment 16 is an overland segment receiving no overland contributions 
itself. Thus columns 1 and 2 of row 16 in IARY have zeros. The lateral 
inflows to segment 39 come from the overland segments 21 and 22 which are 
found in columns 4 and 5 of AUX. Thus columns 1 and 2 of row 39 in IARY 
have the numbers 4 and 5, respectively. 

Columns 3 and 4 of IARY store information about the channel inflows 
for downstream channel segments. If the nth-segment is a channel segment 
receiving channel inflows, columns 4 and 5 of the nth-row in IARY will have 
the column location of these inflow segments in AUX. For example, segment 
12 is an overland flow segment and has no channel inflows, thus columns 3 
and 4 of IARY's row 12 are zeros. Similarly, segment 29 is a channel 
segment with no upstream channel inflows, thus columns 3 and 4 of row 29 in 
IARY are also zeros. Alternatively, segment 41 has one upstream channel 
inflow (segment 40) which is found in column 1 of AUX when computed. 
Therefore, column 3 of row 41 in IARY contains a 1 while column 4 is a 
zero. Segment 39 has two channel inflows (segments 37 and 38). When these 
upstream segments are computed they are entered respectively in columns 2 
and 3 of AUX. Therefore, columns 3 and 4 of row 39 in IARY contain 2 and 
3, respectively. 

Column 5 of IARY contains the column location in AUX of the computed 
segment outflows. For example, the outflow from segment 16 when computed 
is entered in column 5 of AUX. Thus row 16, column 5, of IARY contains a 
5. Similarly, the outflow of segment 28 is found in column 3 of AUX. Thus 
row 28, column 5, of IARY contains a 3. 

The matrices ISEG and IARY are used sequentially to route water and 
sediment through the flow path for each time step. For each segment number 
listed in ISEG, the numbers in the corresponding row of IARY tell the 
program where the associated inflows are stored in Q and GS. When all the 
segments have been worked through, the entire procedure is repeated for the 
next time step. 
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A1.3 FORTRAN VARIABLES NOT INCLUDED IN THE INPUT DATA LIST 


Name 

ACCLN 

ADP 

ALP 

AP 

BEG 

CND 

CONC 

CP 

DFT 

DX 

END 


ERO 

GCD 

GDLAT 

GDUP 

GS 


GTOT 

HYR 


Description 

Gravitational acceleration 

Thickness of detached soil at the start of the 
current time step 

Kinematic-wave friction coefficient in Eq. 26 
Cross-sectional area of flow at the start of the 
current time step 

Distance of space increment used in sediment 
routing to the upstream boundary of the segment 
Canopy cover density 

Volume concentration of individual material 
fraction at carrying capacity 
Volume concentration of individual material 
fraction at the start of the current time step 
Medium size of se . uit fraction 
Space increment used in sediment routing 
Distance travelled by sediment characteristic at 
the end of the current time step, and measured 
from upstream boundary of segment 

Volume rate of soil detachment by raindrop impact 
Ground cover density 

Volume rate of lateral inflow of material fraction 
to segment during current time step 
Volume rate of upstream inflow of material 
fraction to segment 

Volume discharge of individual material fraction. 
This is a five-column matrix used to store the 
sediment inflow and outflow for a segment. This 
matrix has as many rows as time steps in the 
simulation period. 

Total sediment discharge at the catchment outlet 
Flow depth in overland segments, and hydraulic 
radius in channels 


Units 

ft/sec 2 

ft 

ft^/sec 

ft 2 

ft 

ft 2 /ft 2 

ft 3 /ft 3 

ft 3 /ft 3 

ft 

ft 


ft 

ft 3 /ft? sec 
ft 2 /ft 2 

ft 3 /ft. sec 

ft 3 /sec 


ft 3 /sec 
lbs/sec 

ft 








UPON Number of time steps until time of ponding from 
the start of the simulation period 

ITYPE Index defining type of flow in segment. Set 

ITYPE equal to 1 for overland flow, 2 for channel 
flows with negligible infiltration, and 3 for 
channel flows with significant infiltration 

KOUT Column in matrix IARY containing storage location 

of segment outflow computed at the end of the 
current time step. 

KS Number of time steps until the time of formation 

of a characteristic, or shock, wave. 

K! Number of the current time step 

K2 Index identifying a characteristic wave emanating 

from dry ground (K2=0), or from upstream boundary 
(K2>0) 

NO Number of current space increment when routing 

a sediment characteristic 

NSEG Total number of segments in the catchment 

PERM Saturated hydraulic conductivity in/hr 

POR Porosity of bed material ft 3 /ft 3 

Q Water discharge. This is a five-column matrix 

used to store the water inflow and outflow 
computed for a segment. This matrix has as many 


rows as time steps in the simulation period. ft 3 /sec 

QE Water discharge at the end of the current sediment 

routing step. ft 3 /sec 

QINF Rate of infiltration during the current time step in/hr 

QLAT Lateral inflow of water to segment during the 

current time step. ft 3 /ft. sec 

QOUT Water discharge at the catchment outlet ft 3 /sec 

QUP Upstream inflow of water to segment ft 3 /sec 

RAIN Rain precipitation inches 

RNET Net rainfall intensity in/hr 

S Specific gravity of sediment 

TC Critical bed tractive force lbs/ft 2 
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UST Bed shear velocity 

VEL Average velocity of flow 

VINT Total interception per unit area of land 

VS Settling velocity in quiescent water for material 

fraction 

X Distance measured from upstream boundary of segment 

XS Distance travelled by a characteristic or shock 

water wave during the current time step 
XSS Space increment used for sediment routing (Eq. 81) 

ZD Thickness of detached soil for any one material 

fraction 
Bed elevation 


ft/sec 

ft/sec 

inches 

ft/sec 

ft 

ft 

ft 


ZL 


ft 

ft 
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PROGRAM SEDLAB 

C ************** 

c 

DIMENSION TITLE <20 > .SEG<100),1 SEC(100),QOUT<300),GTOT(300>, 
*SLEN<100>,SLOPE(100),CPER(i00 >,EPER<100),IARY(100,5),OMESI300), 
♦GMES<300 ) ,ZZ<3>,C<?1> ,OVA< 100) ,NEED<6) .CANOUOO) ,GCOV< 100 > , 
*HYCND<100 >,SORPTY<100), STORM!5) 

INTEGER SEG 

COMMON /GEN/ UMAX,ITGOM,ITPON,DT,DTS 

COMMON /COV/ CND.GCD,HLR,HIR 

COMMON /RAI/ DR<300),RNET < 300 >,QINF<300) 

COMMON /INT/ VOR,VOG,SRG,VIN,EVP,VINT,RAIN 
COMMON /INF/ PERM,SORP 

COMMON /FLO/ 0<300,S),ALP,EXP,K,ITYPE,SLN,SLP,CPR,EPR, 
*KOUT,QUP<300 > ,(3LAT<300> 

COMMON /SED/ GS<300,5,5),ZL<100>,ZD<100,S),CDUP<300,S), 

♦GDLAT(300,5) ,POR <S) ,CP <5>,DSO<S),PC<S>,VS<5>,ERD<300), 

♦ BEGQOO) ,END!100 > , ADF , SNU , EMV, HOC , SGC , QF., AP ,K1,BRE,APS,ITS,X, 
*XS,XSS,NO,BE T,DX,SUBW,SPGR,ACCLN,GAMA,K 2,K S,GMAX,NSOIL 
DATA ZZ/'I','M','C'/,DEC/'.'/.BLANK/' '/ 

C 

C DATA INPUT 
C 

READ!S,SO 1) TITLE 
C WATERSHED GEOMETRY 

READ!S,502) AREA,NOV,HCH,NCI,NSTRM 

NTO=NOV+NCH 

NSEG=NTO+NCI 

IC=NOV+l 

READ< 5,503) <SEG<1),OVA<I>,SLEN<I),SLOPE<I),CPER<I),EPER<I), 
*1=1,NOV) 

READ(5,504) < SEG<1>,SL EN<I),SlOPE <I),CPER(I>,EPER<I>,I=IC,NSEG> 
C VEGETATIVE COVER AND SOIL CHARACTERISTICS 
READ< 5 , SOS) VOG ,SP.G , VOR ,HLR 
READ!5,S06) NSOIL.SPGR,AIM,ADF,GMAX,GDX 
READ<5,516) <DS0 < I ) , I---1 ,NSOIL ) 

READ <5,517) <PC<I),I = 1,NSOIL> 

READ<5,507) <CANO<I),CCOV<I),HYCND<I), 

*1=1,NSEG) 

C COMPUTATIONAL SEQUENCE 

READ< 5,500) < 1SEC < I) , < IAR Y < I , J ) , J = 1 , S) , I --1 , NSEG) 

C WATER PROPERTIES AND FLOW ROUTING PARAMETERS 
READ<5,509) TEMP,GAMA,SNU ,EXP,C1,C2 
C OUTPUT OPTIONS 

READ < S,510) <NEED<I),1 = 1,6) 

C 

C LISTING OF INPUT DATA 
C 

WRITE< 6,601) TITLF 

IF < NEED <1) EQ 0) GO TO 9S 

WRITE <6,603) AREA , NOV , NC.H , NCI , NSEG , NSTRM 

WRITE<6,604) <SEG<I),OVA<I>,SLEN(I),SLOPE<I),CPER<I),EPER<I), 
*CANO<I),GCOV(I),HYCND <I),1 = 1,NOV> 

WRITE<6,605) <SEG<I),SLEN <I),SLOPE!I),CPER<I),EPER<I),CANO<I), 
*GCOV < I) , HYCND < I ) , I--IC.NSEG) 

WRITE <6,606) VOG , SRC,, VOR , HLR , NSOIL , SPGR , AIM, ADF , GMAX, GDX 

WRITE < 6,616) 

WRITE<6,617) <I,D50<I),PC<I>,I=1,NSOIL) 

WRITE<6,618) 

WRITE<6,607) <ISEG<I),SEG<I),<IARY(I,J>,J=1,5),1=1,NSEG> 
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URITE(6,606) TEMP,GAMA,SNU,EXP,Ci,C2 
95 CONTINUE 
C 

C INVARIANT INFORMATION 
ACCLN=32.2 

SUBW=(SPGR-1.0)*GAMA 
C POROSITY AND SETTLING VELOCITY 
DO 88 1=1,NSOIL 

POR(I)=1.-0.245-0.0864/(0.1*D50<I)>**0.21 
DMM=D50(I> 

CALL SETVEL<DMM,TEMP,SPGR,W) 

VS(I)=W 

D50(I)=D50(I)/304.8 
88 CONTINUE 
C 

TRAIN=0.0 
TVINT=0.0 
C 

C COMPUTATION FOR EACH STORM 
DO 100 N=1,NSTRM 
C 

C STORM DATA INPUT 

READ!S,SI1) STORM,DTM,ITMAX,ITCOM,EVP,VIN 
READ < S,S12) (SORP TY <I),I = 1,NSEG) 

READ(S,513) <DR(I>,1=1,ITMAX) 

READ < S,S14) <QMES(I),1 = 1,ITCOM) 

READ(S,515) (GMES(I),1=1,ITCOM) 

DT=DTM/60.0 
DTS=DTM*60.0 
C 

C LISTING OF STORM DATA 

IF<NEED<1) .GT.0 ) URITE(6,6ii> STORM,VIN,EVP,DTM,ITMAX,ITCOM 
IF(NEED<1).GT.0) URITE<6,612> <SORPTY(I>,1=1,NSEG> 

C 

C POTENTIAL EROSION BY RAINDROP IMPACT 
DO 130 IT=1,ITCOM 
IFUT.GT. ITMAX) DR(IT) = 0.0 
ERO(IT)=AIM*(DR(IT)/4320 0.0)**2.0 
130 CONTINUE 
C 

C ROUTING SEGMENTS ACCORDING TO COMPUTATIONAL SEQUENCE 
URITE(6,621) 

IF(NEED<2).GT.0) WRITE<6,648> 

DO 101 1=1,NSEG 
C 

K=ISEG<I) 

IF(K.LE.NOV) ITYPE=i 

IF< K.GT NOV.AND.K.LE.NTO) ITYPE = 2 

IF < K,GT NTO) ITYPE=3 

SLN=SLEN< K) 

SLP =SLOPE < K) 

CPR=CPER(K) 

EPR=EPER <K) 

CND=CANO < K > 

GCD=GCOV(K) 

PERM=HYCND(K) 

SORP=SORPTY(K) 

EMV=HLR#<1 0-GCD) 

IF < GCD.EQ.0,0) EMV=0.0 
HIR=0.0 
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HGC=HIR*GCD 

5GC=GCD-HGC 

E»RE=< 1 0 -CND) * < 1 0 -CCD ) 

QVF=OUA(K>/<AREA*43560 0) 

C 

C SELECT SPATIAL INCREMENT FOR SEDIMENT ROUTING 
RDX=SLN/GDX 
NDX=INT(RDX) 

IF((RDX-NDX).GT.0.0 ) NDX=NDX+1 
IF < NDX. LT.2) NDX = 2 
DX=SLN/FLOAT(NDX) 

C KINEMATIC FRICTION PARAMETERS 
ALP=C1+C2*SLP**0 5 
C 

C INITIALIZE SURFACE ELEUATION AND DETACHED SOIL STORAGE 
DO 107 J=i,100 
Zl.< J)=0.0 
DO 87 L=i,NSOIL 
ZD(J,L > = 0 0 
87 CONTINUE 

107 CONTINUE 

IF(ITYPE.ER.2) GO TO 108 
C 

C COMPUTING NET RAINFALL RATE 
CALL INTCPN 

IF(ITYPE.EQ i) TRAIN=TRAJN+RAIN*OUF 
IF <ITYPE.EQ.1) TVINT-TVINT+UIN7#OUF 
C 

C COMPUTING PONDING TIME AND POTENTIAL INFILTRATION 
CALL INFLTN 

IF ( UPON . EQ . 0 ) GO TO 151 

108 CONTINUE 
C 

C WATER AND SEDIMENT ROUTING 

C 

C UPSTREAM INI LOU AND LATERAL INFLOU/OUTFLOW 
DO 102 IT=ITPON,ITCOM 
QUP(IT) = 0.0 
QLAT< IT)==0 0 
DO 81 L=1,NSOIL 
GDUP< IT , L > = 0 .0 

81 GDLAT(IT >L)~0.0 

IF<IARY(K,3).EQ.0) GO TO 103 
DO 104 J=-l ,2 

IF(IARY<K,2+J) EQ 0) GO TO 103 
J J = IARY(K,J+2) 

QUP(IT)=QUP(IT)+Q<IT,JJ> 

DO 82 L=i,NSOIL 

82 GDUP(IT,L)= GDUP(IT,L )+GS(IT , J J,L ) 

104 CONTINUE 

103 IF<ITYPE.GT 1) GO TO 505 

QLAT <IT)=QLAT(IT) + (RNET<IT)-QINF(IT))/4320 0.0 

105 IF(IARY(K,1 ) . EQ.0 > GO TO 102 
DO 106 J = 1 ,? 

IF(IARY(K,J).EQ.O) GO TO 102 
JJ=IARY(K,J) 

QLAT ( IT ) - Q L. A T ( I T )+Q( IT ,JJ> 

DO 83 L~1,NSOIL 

83 GDLAT (IT , 1.) - GDLAT < IT , L ) ■* GS < IT , J J , L > 

106 CONTINUE 








102 CONTINUE 

KOUT = I ARY < K,S> 

C 

CALL UROUT 
C 

C LISTING THE FINAL SURFACE LEVEL OF THE SEGMENT 
IF(NEED(2 >.GT.0 > URITE(6,649> K 

IF <NE£D<2> .GT 0) URITE(6,650> (BEG(J),END<J>,ZL(J>,J*1,NO) 

C 

101 CONTINUE 
C 

C RUNOFF AND SEDIMENT DISCHARGE AT THE CATCHMENT OUTLET 
DO 121 IT * 1,1TCOM 
IF<IT.C£ ITPON) GO TO 122 
QOUT ( I T > = 0 0 
GTOT(IT)*0 0 
GO TO 121 

122 QOUT( 1T ) =Q(IT, 1 > 

GTOT <IT > = 0 0 

DO 84 L=l,NSOIL 

84 GTOT <IT)-GTOT(IT)«CS(IT,1,L > *GAMA 
121 CONTINUE 

UR ITE(6,622) 

C 

C COMPUTED YIELDS, PEAKS AND TIME TO PEAKS OF WATER AND SEDIMENT 
TQOUT =0 0 
QCPK=0 0 
TGOUT = 0 0 
GCPK=0 0 

DO 123 IT = ITPON,ITCOM 

IF(QOUT(IT).LT . QCPK) GO TO 124 

QCPK=QOUT(IT) 

IQC=IT 

124 TQOUT = TQOUT + QOUT <IT)#DTS 

IF(GTOT(IT).LT.GCPK) GO TO 125 
GCPK=GTOT(IT) 

IGC=IT 

125 TGOUT =TGOUT+GTOT(IT)*DTS 

123 CONTINUE 

TQPSF=TQOUT/(AREA*43560.0>*12 0 
GO TO 126 
151 CONTINUE 
C 

C SMALL EVENT GENERATING NO RUNOFF 
TQOUT=0.0 
TQPSF=0.0 
QCPK = 0 . 0 
IQC = 0 
TGOUT=0.0 
GCPK=0.0 
IGC = 0 

126 CONTINUE 

IF(NEED<3).EQ.0) GO TO 96 
C 

C WATER BUDGET 

TINF=TRAIN-(TVINT+TQPSF> 

C 

C MEASURED YIELDS, PEAKS AND TIME TO PEAKS OF WATER AND SEDIMENT 
TQMES=0.0 
QMPK=0.0 
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TGMES=0.0 
GMPK-0 0 

DO 127 IT = 1 , ITCOM 
IF<QM£S(IT) LT.QMPK) GO TO 128 
QMPK=QM£S<IT) 

IQM=IT 

128 TQMES~TQMES+QMES< IT)X<DTS 
IF(GMES(IT).LT.GMPK) GO TO 129 
GMPK=GMES(IT) 

IGM=IT 

129 TGMES=TGMES+GMES(IT)*DTS 
127 CONTINUE 

PERCENTAGE DIFFERENCE OF COMPUTED VALUES TO MEASURED VALUES 
PCWY = < TQMES-TQOUT)/TQMES*l00.0 
PCWP=(QMPK~QCPK)/QMPK*100 .0 
PCWT=FLOAT(IQM-IQC)/FLOAT(IQM>*100 . 0 
PCSY=<TGM£S-TGOUT) /TGMF.S* 1 0 0 . 0 
PCSP= < GMPK-GCPK)/GhPK$l0D.0 
PCST = FLOAT<IGM-IGC>/FLOAT(IGM>*100 . 0 

OUTPUT OF WATER BUDGET 

WRITE(6,6S1) TRAIN,TVINT,TINF,TQPSF 

OUTPUT OF YIELDS, PEAKS AND TIMF TO PEAKS OF WATER AND SEDIMENT 

WRITE(6,652) TQMES,TQOUT,PCWY,QMPK.QCPK,PCWP,IQM,IQC,PCWT,TGMES, 
♦TGOUT,PCSY,GMPK,GCPK,PCSP,IGM,IGC,PCST 
CONTINUE 

IF(NEED(4).EQ.0) GO TO 115 

OUTPUT OF RAINFALL INTENSITIES, AND MEASURED AND COMPUTED HYDROGRAPHS 
WRITE < 6,656) 

IP=ITPON-l 

WRITE (6,657) < J , DR < J ) ,RNF.T < J ) , QMES (J ) ,QOUT(J> ,GMES(J) ,GTOT( J) , 

*J=1,IP) 

WRITE<6,658) <J,DR(J),RNET(J),QINF(J),QMES(J),QOUT(J),GMES(J), 
*GTOT<J),J=ITPON,ITCOM) 

115 CONTINUE 
C 

C OUTPUT OF HEITOGRAPH, MEASURED AND COMPUTED HYDROGRAPH 
IF(NEED < 5 > .EQ.0) GO TO 97 
UR ITE(6,66 0 ) 

QMAX= QMPK 

IF(QCPK.GT.QMPK) QMAX=QCPK 
DO 141 J=1,71 

141 G(J)=BLANK 

DO 142 J"1,ITCOM , 2 
G(1 ) = DEC 
IR=DR(J)*7.0 
IRR=71 -IR 
DO 143 J1=IRR,71 

143 G(J1)=ZZ<1> 

IM=QMES<J)/QMAX*70.0 
G <IM > =ZZ < 2) 

IC=QOUT<J)/QMAX*70.0 
G<IC)-ZZ(3 > 

WRITE < 6,661)J,G 
DO 144 Jl-1,71 

144 G < Ji)=BLANK 

142 CONTINUE 









97 CONTINUE 
C 

C OUTPUT OF HEITOGRAPH, MEASURED AND COMPUTED SEDIMENTGRAPH 
IF<NEED(6).EQ.0) GO TO 98 
WRITE<6,662> 

GMAX=GMPK 

IF (GCPK . GT . GMPX ) GMAX«=GCPK 
DO 161 J=1,71 

161 G<J)=BLANK 

DO 162 J=1,ITC0M,2 
G(1 )*DEC 
IR*DR<J>*7.0 
IRR- 71-IR 
DO 163 Ji=IRR,71 

163 G<J1)=ZZ<1) 

IM=GMES<J)/GMAX*70.0 
G(IM)=ZZ < 2) 

IC=GTOT<J)/GMAX*70.0 
G(IC)=ZZ(3> 

WRITE(6,663 >J,G 
DO 164 Jl=l,71 

164 G<Ji)=BLANK 

162 CONTINUE 

98 CONTINUE 
WRITE(6,65?) 

100 CONTINUE 
C 

SOI FORMAT <20A4) 

SU2 FORMAT'FiO.4,414) 

503 FORMAT(I4,F12.3,4F10.4) 

504 FORMAT<I4,4F10.4) 

505 FORMAT(5F10.4) 

506 FORMAT(I4,F10,4,2F10 8,2F10.4) 

516 FORMAT(10F7.3 > 

517 FORMAT(10F7.3) 

507 FORMAT(3F10.4) 

508 FORMAT <614) 

509 FORMAT(2F10.4,F10.8,3F10.4> 

510 FORMAT(614) 

511 FORMAT(SA4,F10.4,214,2F10.4) 

512 FORMAT < 8F7.4) 

513 FORMAT(12F6.3) 

514 FORMAT(12F6.2) 

515 FORMAT <12F6.2) 

601 FORMAT(1H1///X,20A4///) 

603 FORMAT(30X,'CATCHMENT GEOMETRY'//4IX,' DRAINAGE AREA = ', 

♦F7.1,' ACRES '/38X,'OVERLAND SEGMENTS = ',IS/10X,'CHANNEL ', 
♦'SEGMENTS WITH NEGLIGIBLE INFILTRATION = ',I5/8XCHANNEL ', 
♦'SEGMENTS WITH CONSIDERABLE INFILTRATION = ',15/41X,'TOTAL ', 
♦'SEGMENTS = ',IS//30X,'NUMBER OF COMPUTED STORMS = ',15//// 

♦24X,'PHYSICAL CHARACTERISTICS OF SEGMENTS'// X,'- 

♦---' 

♦/36X,'RELATION BETWEEN CANDPY GROUND SATURATED'/X,'SEGMENT OVE 
♦RLAND LENGTH SLOPE WETTED PERIMETER COVER COVER HYDRAULIC 
♦'/X,'NUMBER AREA',21X,'AND FLOW AREA DENS. DENSITY CONDU 

♦CT.'/X,'4 TYPE <SQ FT) <FT> <FT/FT) COEFF. EXP.',20X, 

♦ '( INCH/HR >'/X, '-' 

♦ ,'-'//) 

604 FORMAT<X,13,' OV',X,F10.1,X,F7.1, 2X,F7.5,2X,F5.1,3X,F5.3,4X,F5.3, 
♦3X,F5.3,4X,F5.3) 











605 FORMAT<X,13,' CH',13X,F7 1,X,F7.S,2X,FS.1,3X,FS.3,4X,FS.3,3X, 
tF5.3^X,FS 3) 

606 FORMAT<////27X,'GROUND COVER CHARACTER 1ST ICS'//6X, 

♦ 'INTERCEPTION STORAGE CAPACITY FOR GROUND COVER ■= ', 

♦F9 3,' INCHES',/6X,'RATIO OF EVAPORATING SURFACE', 

♦ ' TO PROJECTED AREA = ',F9.3/7X,'RAT 10 BETWEEN INTERCEPTION'. 

♦ ' STORAGE CAPACITIESV19X , 'OF CANOPY COVER AND GROUND COVER', 

♦ ' - ',F9.3/10X,'AVERAGE HEIGHT OF GROUND COVER IN CHANNELS', 

♦ ' = ' ,F9.3, ' FEET' 

♦ Z///18X,'SOIL CHARACTERISTICS AND DETACHMENT PARAMETERS'// 
♦27X,'NUMBER OF SEDIMENT SIZE FRACTIONS = ',15/ 

♦32X,'SPECIFIC GRAVITY OF SEDIMENT - ',FB.2/liX, 

*'COEFFICIENT OF SOIL DETACHMENT BY RAINDROP IMPACT = ', 

♦E10 3/17X,'COEFFICIENT OF SOIL EROSION BY SURFACE FLOW = 

♦E10.3/22X,'MAXIMUN PENETRATION DEPTH UF RAINDROPS = ',F9.3, 

♦ ' FEET'/20X,'SPACE INCREMENT USED IN SEDIMENT ROUTING', 

♦ ' = ',F8.2,'FEET'//> 

616 FORMAT <1SX,'- 

♦/16X,'SERIAL NO. OF MEDIUM SIZE OF PERCENTAGE'/I6X, 

♦ 'SED FRACTION SEDIMENT FRACTION OF SEDIMENT'/37X, 

♦' <MM)' ,13X, 'FRACTI0N'/15X, '- 

♦ ' - - . /r> 

617 FORMAT(19X,I4,l3X,F7.3,12X,FS.i> 

618 FORMAT <////SX, '- - 

♦ ,'-'/9X, 'ISEG' ,30X, ' IARYV27X,'THIS ARRAY ' 

♦'CONTAINS THE STORAGE LOCATIONS'/SX,'COMPUTATIONAL',9X,'OF ', 
♦'COMPUTED INFLOWS AND OUTFLOWS IN THE'/8X,'SEQUENCE1IX,'WATER' 

♦ '<(■>) AND SEDIMENT < GS> DISCHARGE MATRICES'/ 

♦ 22X , '-'/ 

♦ 11X,'♦',1 OX,'SEGMENT t LATERAL INFLOWS UPSTREAM INFLOWS ', 

♦ 'OUTFLOU'/SX,'- - - ', 

♦ '- -'//) 

607 FORMAT<1 OX,13,11X,13,9X,11,6X,11,9X,11,8X,11,BX,11> 

608 F0RMAT<////18X,'WATER PROPERTIES AND KINEMATIC FRICTION 
♦'PARAMETERS'//31X,'WATER TEMPERATURE =' 

♦ ,F8.2 , ' I)EG CENT '/24X,'SPECIFIC WEIGHT OF WATER *' 

♦,F8.2,' LBS/CU.KT '/20X,'KINEMATIC VISCOSITY OF WATER = ', 

♦ HX , El 0.3,' SQ.FT./SEC.'/13X, 

♦'EXPONENT IN KINEMATIC WAVE EQUATION = ',F7 2/19X, 

♦'KINEMATIC FRICTION PARAMETERS , C1=',F6.2,', C2 =',F6.2//> 

611 FORMAT(//36X,'STORM DATA'//37X,'STORM DATE = ',5A4/16X, 

♦ 'ANTECEDENT INTF RCEPTION STORAGE = ',F10.4,' INCHES'/26X, 

♦'MEAN EVAPORATION RATE = ',F10.4,' INCHES PER HOUR'// 

♦ 2SX,'TI ME STEP DURING STORM *',F8.1,' MINUTES'/ 

♦8X,'NUMBER OF TIME STEPS IN RAINFALL PERIOD = ', 

♦I5/6X,'NUMBER OF TIME STEPS IN SIMULATION PERIOD = ', 

♦I5///21X,'SORPTIVITY FOR EACH CATCHMENT SEGMENT', 

♦/28X,'(SQUARE INCHES PER HOUR)'/) 

612 FORMAT(X,1 OF 8.4 > 

621 FORMAT(////20X,'***###***#*****##*#*#***#*********####*##**'/ 

♦20X,'* START OF WATER AND SEDIMENT ROUTING *' 

♦ /20X, '*****#********#*###*##**#*#**)***##****#**#*') 

648 FORMAT<////?X, '-' , 

% '-'/9X 'DISTANCE' 

♦' FROM U/S END DISTANCE FROM U/S END CHANGE OF BED ELEV.'/ 
♦8X,'0F SEGMENT TO START OF OF SEGMENT TO END OF ', 

♦'FOR THE SPACE'/iOX,'SPACE INCREMENT (FT)',2X,'SPACE 

♦ 'INCREMENT (FT)',SX,'INCREMFNT<FT)' /9X , '- 

♦ ' - - //> 

649 F0RMAT<36X,'FOR SEGMENT ',14) 





















650 FORMAT (1SX.F7.1 , 17X,F7 1 , 1SX,F6.3> 

622 FORMAT <////15X,'a*******************************************', 

♦'*»*#***'/15X, '* END OF WATER AND SEDIMENT ROUTING', 

♦ ' *'/15X,'*******************»#********«', 

S'**********#**##*****#'//) 

651 F0RMAT(//35X,'WATER BUDGET'//32X,'PRECIPITATION = ',F8.2, 

♦ ' INCHES'/13X,'LOSS DUE TO INTERCEPTION STORAGE = ',F8.2, 

♦ ' INCHESV21X,'LOSS DUE TO INFILTRATION - ',F8.2,' INCHES'/ 

♦34X,'WATER YIELD = ',F8.2,' INCHES'///) 

652 FORMAT<//33X, 'SUMMARY OF RESULTS'/X,'-', 

♦ ' - ' , 

♦ '-'/8X,'ITEM ',1SX,'MEASURED',7X,'COMPUTED',9X,'UNITS',8X, 

♦ 'PER CENT ' /73X , 'DIFFER . '/X, '-' ,2X, '-' , 

♦ '-' ,2X,'- - -'// 

♦ 6X,'WATER YIELD',7X,F12.1,3X,F12.1,6X,'CUBIC FEET',6X,F6 1,/X, 

♦'WATER DISCHARGE PEAK' ,7X,F8 1,7X,F8.1,4X,'CUBIC FEET/SEC', 
♦4X,F6.i/2X,'TIME TO WATER P F.AK ' , 9X , 15,1 OX, IS , 9X, ' TIME STEPS', 
♦5X,F6.1/4X,'SEDIMENT YIELD',6X,F12. 1,3X,F12.1,8X,'POUNDS',8X, 
♦F6.1/2X,'SEDIM. DISCHARGE PEAK',5X,FB.1,7X,F8.1,5X,'POUNDS ', 
♦'PER SEC', F9.1/X,'TIME TO SEDIMENT PEAK',7X,15,10X,IS, 

♦9X , 'TIME STEPS',SX,F6.l//> 

656 FORMAT (///X, '-', 

♦'-'/X, 'TIME' ,3X ,' PRECIP ' 

♦,SX,'NET',6X,'POTENTIAL',9X,'RUNOFF',8X,'SEDIMENT DISCHARGE'/X, 

♦ 'STEP',11X,'RAINFALL INFILTR AT ION' , 2X , '-' 

♦ ,2X, '-'/40X, 'MEASURED COMPUTED MEASURED ', 

♦'COMPUTED'/3X,'#',SX,'IN/HR',3X,'IN./HR.',SX,'IN./HR ' ,7X, 

♦'CFS',7X,'CFS',5X,'LBS/SEC',3X,'LBS/SEC'/X,'- - ', 

♦ - - - - - -' 

♦ /) 

657 FORMAT<X,IA,3X,F6 3,3X,F6.3,8X,'-',BX,F7.1,3X,F7.1, 

♦3X,F7 1,3X,F7.1) 

658 FORMAT<X,I4,3X,F6.3,3X,F6.3,6X,F6.3,5X,F7.1,3X,F7.1, 

♦3X,F7.1,3X,F7.1) 

659 F0RMAT<//8<'**********')) 

660 F0RMAT<///4X,7('♦****♦****'),'*'/4X,'* PLOTS OF HYETOGRAPH< I) ■, ' 

♦ ' COMPUTED<C) AND MEASURED<M> HYDROGRAPHS *' /4X, 


♦7('****»**Ht**'),'*'//14X,'(Q/QMAX)*i00 ->',19X,'<- PRECIP', 

♦'(INCHES/HOUR)'/2X,'TIME'/2X,'STEP',37X,'5 4 3 

♦' 1 0'/3X,'#',4X,'0 10 20 30 40 

♦ ' 60 70 80 90 10 0' /2X, '- I',10<'-1')) 

661 FORMAT(2X,I4,2X,71A1) 

662 F0RMAK///4X ,7< '*****#****'),'#'/4X, '* PLOTS OF HYETOGRAPH< I) ; ' 

♦ ' COMPUTED<C) AND MEASURED(M> SEDIGRAPHS *' /4X, 

♦ 7< '*****#**##'), ')K'//14X, ' <G/GMAX)!K100->',19X,'<-PRECIP', 

♦'(INCHES/HOUR)'/2X,'TIME'/2X,'STEP',37X,'5 4 3 2', 

♦ ' 1 0'/3X,'♦',4X,'0 10 20 30 40 50' 

♦ ' 60 70 80 90 100' /2X,'- I',10<'- 1')> 


663 FORMAT <2X,I4,2X,71A1) 

END FILE 6 

STOP 

END 

SUBROUTINE INTCPN 

C THIS SUBROUTINE COMPUTES INTERCEPTION LOSSES AND NET RAINFALL 
C INTENSITY 

COMMON /GEN/ ITMAX,ITCOM,ITPON,DT,DTS 

COMMON /COV/ CND,GCD,HLR,HIR 

COMMON /RAI/ DR(300),RNET<300),QINF<300> 

COMMON /INT/ VOR,UOG,SRG,VIN,EVP,UINT,RAIN 
C 
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C INITIALIZING.CUMMULATIVE VALUES AND STORAGE CAPACITIES 
TR1=0.0 
TR2=0.0 
TR3=0.0 
TR4=0.0 
TRC=0.0 
TRG=0.0 
VOC=VOR*VOG 
SRC=VOR*SRG 
CCAP = < I . 0-VIN)*VOC 
GCAP=(i.0-VIN)*VOG 

COMPUTING NET RAINFALL FOR EACH TIME INTERVAL 
DO 201 IT=1 , ITCOM 

RATE OF INTERCEPTION ON CANOPY AND THROUGHFALL FROM CANOPY 
TR1=TR1+DR(IT)*DT 
IF < TRi.LE.CCAP) DRC = DR (IT > 

IF< TR1 . GT . CCAP > DRC=--EVP*SRC 

IF(TRi GT.CCAP AND TRC.LT.CCAP ) DRC=(CCAP-TRC)/DT 
TRC = TRC-»DRC*DT 
RTHO=DR<IT>-CND#DRC 

RATE OF INTERCEPTION ON GROUND COVER AND NET RAINFALL RATE 
TR2=TR2+RTH0*DT 
IF < TR2.LE.GCAP > DRG = RTHD 
IF<TR2.GT.GCAP> DRG=EVP*SRG 

IF(TR2 .GT CCAP.AND.TRG.LT.GCAP > DRG~(GCAP-TRG)/DT 
TRG--TRG+DRG*DT 
RNET(IT)=RTHO-GOD*DRG 
TR3=TR3+RNET <IT)*D1 
01 CONTINUE 

TOTAL RAINFALL AND TOTAL INTERCEPTION STORAGE 
RAIN-TR1 
VINT=TRi-TR3 
RETURN 
END 

SUBROUTINE INFLTN 

THIS SUBROUTINE COMPUTES THE PONDING TIME AND THE DECAY OF 
POTENTIAL INFILTRATION FROM THE TIME OF PONDING 

COMMON /GEN/ ITMAX,ITCOM,ITPON,DT,DTS 
COMMON /RAI/ DR(300>,RNET<300),QINF(300> 

COMMON /INF/ PERM,SORP 

PONDING TIME 
SL=0.0 
ITPON=0 

DO 251 IT = 1 ,ITCOM 
IF <ITPON.GT.0 > GO TO 252 
SL=SL+RNET(IT)#DT 
IF(RNET(IT) LE.PERM) GO TO 251 
SR-SORP/(RNET(IT)-PERM) 

IF(SL.LT.SR) GO TO 251 
ITPON-IT 
RP=RNET(IT) 

FP“SL 
AIN=RP 
C=RP-PERM 
C1=C/PERM 
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C2=C/RP 
GO TO 253 
C 

C POTENTIAL INFILTRATION, THE DECAY EQUATION IS SOLVED BY NEWTON'S 
C METHOD 

252 SATI=PERM#FLOAT <IT-ITPON)*DT 
ITR = 0 

2S A IF<AIN.LE.PERM) GO TO 253 

C3=C2*AIN/(AIN-PERM) 

FF=SATI-FP*< <RP-AIN>/<AIN-PERM>-Cl*ALOG<C3)) 

DFF=FP*C*<1 0-<AIN-PERM)/AIN)/CAIN-PERM)**2 0 

AINN=AIN-FF/DFF 

TOL=ABS< <AINN-AIN)/AINN> 

AIN=AINN 

ITR=ITR+1 

IFOTR.GT. 100) GO TO 25S 
IF(TOL.GT.0.Oi) GO TO 254 

253 IF<AIN.LE.PERM) AIN=PERM 
QINF<IT)=AIN 

2S1 CONTINUE 
GO TO 256 

255 WRITE < 6,261) 

STOP 

256 RETURN 

261 FORMAT<20X,'********** ITERATION EXCEEDS 100 > INFILTRATION', 
♦ ' EQUATION **********') 

END 

C 

SUBROUTINE UROIIT 
C 

C THIS SUBROUTINE ROUTES WATER FROM OVERLAND AND CHANNEL UNITS 
C THE ROUTING PROCEDURE IS BASED ON THE CHARACTERISTIC SOLUTION 
C OF THE KINEMATIC WAVE APPROXIMATION 

COMMON /GEN/ ITMAX,ITCOM,I TPON,DT,DTS 

COMMON /COV/ CND,GCD,MLR,HIR 

COMMON /RAI/ DR(300),RNET(300),QINF(300) 

COMMON /FLO/ Q<300,5),Al P,EXP,K,I TYPE,SLN,SLP,CPR , EPR , 

♦KOUT,QUP <300),QLAT<300) 

COMMON /SED/ GS(300,5,5) ,ZL(100 ) ,ZD< 10 0 ,5) ,GI)UP (30 0 ,5 ) , 

♦GDLAT(300,5),POR(5),CP<5), D50(5),PCC5),VS(5),ERO(300) , 

♦BEG(100),END(100),ADF,SNU,EMV,HGC,SGC,QE,AP,K1,BRE,APS,ITS,X, 
♦XS,XSS,N0,BET,DX,SUBW,SPCR,ACCLN,GAMA,K2,KS,GMAX,NS0IL 
C 

EXP 1=£XP + i.0 
BET=1.0/EXP 
EXM1=EXP-1.0 
TERM=EXP*ALP*DTS 
C 

Kl=ITPON 
K2"0 

301 IF(K1.GT.ITCOM) GO TO 391 
KS=K 1 
ITS=K1 
NO=0 

QB-QUP < K1) 

AB-(QB/ALP >##BET 

QP = QB 

AP=AB 

APS=AP 

DO 381 L=1,NSOIL 







CP<L>=0 O' 

381 IF(QP GT.O.O) CP<L)~GDUP(K1,L>/QP 
QL=QLAT(K1) 

IF (ITYPE . NE . 3) GO TO 302 
APF=AP+QL*DTS 

IF < APF . GT .0.0) PERIM=CPR*APF**EPR 
IF(APF.LE.0.0) PERIM=0.0 
QL=QL-PERIM*QINF<K11/43200.0 
302 IF<QL.LE.0.0.AND.QB.EQ.0.0) GO TO 303 
X=0 .0 
XSS-0.0 

IF < K2.GT.0 > GO TO 304 
AP=0 .0 
QP=0.0 

DO 382 L=1,NSOIL 

382 CP < L) = 0.0 
GO TO 30S 

304 IF(QB.EQ.0.0) GO TO 305 
C TESTING FOR SHOCK FORMATION 
QA=0.0 

IF(Ki.GT.l) QA=QUP(K1-1> 

AA= < QA/ALP)**BET 
IF <AA.GE.AB) GO TO 30S 
SMOCK-AB-AA 
ADD=QL*DTS*0 5 

IF <K1.GT.1) ADD-ADD+QLAT(Ki-i)*DTS*0.5 
IF(SHOCK.LE.ADD) GO TO 305 
C 

C SOLUTION WITH SHOCK 
311 IF<K1.GT.ITCOM) GO TO 312 
QL-QLAT (K 1 ) 

IF <ITYPE.NT 3) CO TO 313 
APF=AP+QL*D1S 

IF<APF GT 0.0) PERIM=CPR#APF**EPR 
IF< APF . LE . 0 . 0 > PERIM==0 . 0 
QL=QL~PERIM*QINF(Ki)/43200.0 

313 ABF=AB+QL*DTS 

IF(ABF.IE 0.0) ARF = 0 . 0 

IF(ABF EQ.0 0) GO TO 305 

AAF=AA + QL*I)TS 

IF < AAF.LE.0.0) AAF = 0 . 0 

QBF=ALP*ABF**EXP 

QAF=ALP*AAF**EXP 

IF(QL.EQ.0 0) GO TO 314 

DEN=EXP1*< AB-AA)*QL 

PROD=ALP/I)EN 

XS=PROD*(ABF**EXPl-AAF**EXPl-AB**EXPi+AA**EXPl) 
GO TO 322 

314 XS=<QB-QA)*DTS/<AB-AA) 

322 AB=ABF 

AA=AAF 

QB-QBF 

QA-QAF 

X=X+XS 

XSS=XSS+XS 

IF (X. GE . SLN > GO TO 32.3 
C 

C SEDIMENT ROUTING 

QE = ALP*C ( (QB/ALP )**BET+(QA/ALP >#*BET>/2. 0>**F.XP 
IF(XSS.GT.DX) CALL SROUT 
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QP = QB 

AP=<AA+A&>/2.0 
K1=K1+1 
GO TO 311 

323 QB = QF<-QL*<X-SLN) 

QA=QA-QL*(X-SLN> 

QC-ALP * << < QB/ALP)**BET+ <QA/ALP >**BET>/2.0 >**EXP 
IF(K1.EQ KC) QC=(GC+Q(KC,KOUT>)/2.0 
Q(K1,KOUT)=QC 
C 

QE=Q < K1, KOUT) 

CALL SROUT 
DO 390 L=1,NSOIL 
390 GS(K1,KOUT,L)=CP<L)*QE 
GO TO 324 
C 

C FLOW CEASES 
C 

303 IF(K2 EQ.O) GO TO 325 
IF<Ki. GT.KC) GO TO 326 
GO TO 327 

326 K2= 0 

325 Q(K1,KOUT)~0.0 

DO 392 L=1,NSOIL 

392 GS (K1, K CUT,L)=0.0 

327 K1=KS+1 
GO TO 301 

C 

C SOLUTION WITHOUT SHOCK 
30S IF(K1.GT.ITCOM) GO TO 312 
QL = QLATiK1) 

IFdTYPE.NF. 3) CO TO 331 
APF=AP+QL*DTS 

IF (APE' GT 0 0) PERIMd:PR*ArE**EPR 
IF(APF tE 0.0) PER IH-0 0 
QL=QL-PERIM*QINF<K1 )/43200.0 

331 AC = AP +fj)L #DTS 

IF < K 1 FQ.KS AND K2 GT 0> AC=AP+QL*DTS/2.0 

IF<AC.IC 0.0) AC-0 0 

QC=ALP#AC**EXP 

IFIQL.EQ 0.0) CO TO 332 

XS= < QC-QP)/QL 

GO TO 333 

332 XS=TERM*AP**EXM1 

333 X=X+XS 
XSS-XSS+XS 

IF(X CF.5LN AND K2 EQ.O) GO TO 351 
IF(X.GE.SLN) GO TO 334 
IF(K2 CT 0) GO TO 360 
C INITIAL RISING PART OF HYDROGRAPH 
Q<K1,KOUT > = (QC + QP)/2.0 
QE=Q(K1 .KOUT) 

CAl L SROUT 
DO 393 t =1,N5O T L 

393 G S(K1 ,KOU1 ,L)=CP<L)*QE 
GO TO 362 

360 CONTINUE 
C 

C SEDfPFNT ROOTING 
QE = QC 
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IFfXSS.GT.DX> CALL SROUT 
QP=QC 
AP =AC 
Kl=Ki+i 
GO TO 305 
334 Q(K1, KOUT > =QC-QL#<X-SLN) 

C 

Q£=Q<K1,KOUT> 

CALL SROUT 
DO 394 L = 1,NSOIL 

394 GS<Ki,KOUT,L>=CP(L)*QE 
C 

C INTERPOLATION OF FLOW AND SEDIMENT DISCHARGES AT THE SEGMENT OUTLET 
324 K3=Ki-KC 

IF<K3.LE.1) GO TO 3Si 
QD=Q(K1,KOUT)~Q(KC,K OUT) 

QAD=QD/r LClAT (K3> 

K3M=K3-i 
DO 352 J = 1,K3M 

Q(KC+J,K OUT) = Q <KC,KOUT )+QAD*FLOAT<J) 

352 CONTINUE 

DO 395 L~ 1,NSOIL 

GSD=GS<K 1 ,KOUT,L)-GS(KC,KOUT,L > 

GSAD=GSD/F1.0AT<K3> 

DO 396 J=1,K3M 

396 GS(KC+J,KGUT,L)=GS(KC,KOUT,L)+GSAD4FL0AT(J> 

395 CONTINUE 

351 KC = K1 

1F(K2.EQ 0) KC~K1-1 
IF<KC.LT.1) KC=1 
K1>=KS+1 

IF(K2.EQ.0) K i~KS 
K2=KS 
GO TO 301 
C 

C EXTRAPOLATION OF OUTFLOWS AT END OF LAST TIME STEP 
312 IFtKC.EQ TTCOM) GO TO 391 
K3~ITCOM-K C 

QD=Q(KC,KOUT >-Q<KC-1 ; KOUT> 

QAD-'QD/FLOAT < K 3 ) 

DO 353 J = 1 ,K3 

Q(KC+J,KOUT)=Q < KC, KOUT)+QAD#FLOAT(J) 

353 CONTINUE 

DO 397 L=l,NSOIL 

GSD=GS < K C,K OUT,L)-GS(KC-i, KOUT, L) 

GSAD=GSD/F'LOAT < K3) 

DO 390 J = 1 ,K 3 

398 GS<KC+J,KOUT,L)=GS<KC,KOUT,L>+GSAD*FLOAT<J) 

397 CONTINUE 
391 RETURN 

END 

C 

SUBROUTINE SROUT 

C THIS SUBROUTINE ROUTES ALL SEDIMENT FRACTIONS 
DIMENSION ADP<5),GLAT<5>,CONC(S> 

COMMON /GEN/ ITMAX,ITCOM,ITPON , DT,DTS 

COMMON /COV/ CND , GCI), HLR , HIR 

COMMON /RAT/ DR<300 > ,RNET<300>,QINF<300> 

COMMON /FLO/ Q<300,5),ALP,EXP,K,ITYPE,SLN,SLP,CPR,EPR, 

*KOUT,QUP<300>,QLAT <300 ) 
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COMMON /SED/ GS(300 ,5,5) ,ZL< 100 ) ,ZD< 100,5) ,GDUP < 300,5 ) , 
*GDLAT<300,S>,POR(S),CP<S>,D50(5),PC<5),VS<S),ERO<300>, 
♦BEGtlOO),END(100),ADF,SNU,EMV,HGC,SGC,QE,AP,K1,PRE,APS,ITS,X, 
$XS,XSS,NO,BET,DX,SUBW,SPGR,ACCLN,GAMA, K2,K S,GNAX,NSOIL 
C 

C SOIL SURFACE ELEVATION AND DETACHED SOIL VOLUME AT THE START OF 
C THE CURRENT TIME STEP 

IF<K2.EQ.O OR K2.EQ.KS) GO TO 401 
DO 402 J=l,100 

IF(X.GE.BEG< J).AND.X.LE.END(J)> GO TO 403 
GO TO 402 
403 Z=ZL < J) 

DO 490 L=l,NSOIL 

490 ADP(L)= ZD < J,L) 

GO TO 480 

402 CONTINUE 
GO TO 480 

401 IF(K1.EQ.KS > Z = 0.0 

IF <K1.GT.KS) Z=ZL(NO) 

DO 496 L=1,NSOIL 
IF(Kl.EQ.KS) ADP < L ) = 0.0 
496 IF(K1.GT.KS) ADP(L>=ZD(NO,L) 

480 NO=NO+l 

DEG < NO)=X-XSS 
END(NO)=X 

C HYDRAULIC PARAMETERS 

IF(QE.LT.l.OE-S) GO TO 404 
AE= < QE/ALP)*#BET 
VE1.=QE/AE 
UEP=CPR*AE**£PR 
HYR=AE/WEP 
RN=QE/<SNU*WEP ) 

C AVERAGE LATERAL INFLOW AND POTENTIAL EROSION RATES DURING THE 
C CURRENT TIME STEP 
K4=Ki-ITS+l 
DO 491 L--1,NSOIL 
TGLAT=0.0 
DO 4 OS J = ITS , K1 
4 OS TGLAT = TGLAT + GDLAT<J,L> 

491 GLAT <L)=TGLAT/FLOAT < K4) 

TFRO=0 0 

DO 406 J=ITS,K1 
406 TERO”TERO+ERO< J) 

AERO-TERO/FLOAT <K4) 

DTSS=DTS*FLOAT < K 4) 

IF(X GT.SLN) DTSS=DTS/(1.+(X-SLN>/<SLN-X+XS>> 

DIB=BRE*AERO 

C COMPUTE EFFECTIVE TRACTIVE FORCE IN VEGETATED REACHES 
IF<ITYPE EO.1 OR.EMV.Eq.0.0) GO TO 408 
IF(Z GT.EMV) GO TO 409 
RATIO=l.0-Z/EMV 
IF(RATIO.GT.1.0) RATIO=i.0 
EHT=HLR*RATIO 
E5C=SGC*RATI0 
EGC=i.0-HGC-ESC 

IF < HYR.GT.EHT) EGC = 1.0-HGC-ESC*EHT/HYR 
GO TO 410 
409 EGC=1.0-HOC 

GO TO 410 
408 EGC=i.0-GCD 
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410 TAO=GAMA*HYR*SLP*EGC 
UST = SQR T C TACI4ACCLN/GAMA > 

DO 499 L=i ,NSOIL 

C SEDIMENT CARRYING CAPACITY 
DFT=D50 <L) 

DMM=DFT*304 8 
W=VS<L> 

IF{ITYPE.EG).1)CALL YALIN(UST,DFT,SNU,SUBW,SPGR,ACCLN,VEL,HYR, 
♦CEE) 

IFIITYPE.GT i.AND.DMM.GE.O 1) CALL YANG(DFT,UST,SNU,VEL,SLP, 
tU,SPGR,CEE) 

IF(ITYPE.GT.1.AND.DhM. LT . 0 . 1) CALL LAURSN(DFT,UST,SPGR ,VEL, 
♦GAMA,SUBW,SNU,DFT,HYR,ACCLN,W,CEE) 

CONC(L)=CEL*PC(L>/100 . 0 

C POTENTIAL SEDIMENT EROSION(+VE>/DEPOSITIONC-VE> 

GDPL=<AE*CONC<L)-APS*CP <L>>/DTSS-GLAT(L> 

BT=W*DTSS/HYR 

IF(GDPL.GT.0.0 >CO TO 407 

IF< BT.LT . 1.0)GDPL=BT*GDPL 

C DEPTH OK DETACHED SOIL AT THE START OF THE CURRENT TIME STEP 
407 ADC=ADP(L ) 

C EROSION BY RAINDROP IMPACT 
RATIO=l.O-ADC/GMAX 
IF (RATIO . l.E .0.0) RATIO=0 . 0 

IF(RN LE 900.0) ADC=ADC+PC<L)/i00.0*DIB#RAT104DTSS 
RAD=ADC-CDPL/WEP*DTSS 
C EROSION BY SURFACE RUNOFF 

IF(RAD.GE.0.0 > GO TO 411 

DGD=~ADF*RAB 

GD-UEP# < ADC+DGD)/DTSS 

AS=APS*CP<L) + (GD+GLAT<L> >*DTSS 

CGNC(L)=AS/AE 

2D(NO, L ) - 0 0 

DZ=-DGD 

GO TO 413 

C DETACHED SOIL DEPTH AND SOIL SURFACE ELEVATION AT THE END OF THE 
C CURRENT TIME STEP 

411 ZD < NO,L)-RAD 

DZ= (R AD-ADP CD) /F.GC 
413 ZL(NO)-Z + I)7/POR (L) 

412 CP(L)=CONC(L) 

499 CONTINUE 

GO TO 497 
C RUNOFF CEASED 
404 AE=0.0 

QE-0 0 

DO 498 L-'l , NSOIL 
498 CP < L)=0 . 0 

497 APS-AE 

ITS=K1+1 
XSS=0.0 
RETURN 
END 
C 

SUBROUTINE YALIN(UST,DFT,SNU,SUBW,S,G,VEL,HYR,CEE) 

C 

C THIS SUBROUTINE COMPUTES THE TRANSPORT RATE OF NON COHESIVE 
C SEDIMENTS USING THE BEDLOAD FORMULA DEVELOPED BY YALIN<1963> 

C 

CALL SHIELD<UST,DFT,SNU,SUBW,TC) 
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Y=UST*UST/<<S-i.)*G*DFT) 

YC=TC/(SUBU*DFT) 

IF < Y-YC > 4,4,5 

4 CEE=0.0 
GO TO 6 

5 SEGMA=Y/YC-1.0 

A=2.45*SQRT(YC)/S**0 . 4 

X1 = 0.63S*DFT#UST*S*SEGMA 

X2=1.-ALOG(i.+A*SEGMA)/< A*SEGMA) 

CEE=X1*X2/(VEL*HYR> 

6 CONTINUE 
RETURN 
END 

C 

SUBROUTINE YANG < DFT , UST , SNU, V, SL.P , W, S, CEE ) 

C 

C THIS ROUTINE COMPUTES THE TRANSPORT RATE OF NONCOHESIVE SEDIMENTS 
C USING THE TOTAL LOAD FORMULA DEVELOPED BY YANG<1973) 

D=DFT 

A=UST*D/SNU 

IF < A.GE.70.0) GO TO 7 

VCU=2.5/< ALOG10 < A)-0.06> + 0.66 

GO TO 8 

VCW-2.05 

ESP=V*SLP/W-VCW*SLP 
IF (ESP) 9,9,10 
CEE = 0.0 
GO TO 11 

0 FI=5.435 -0 286*AL0G10 (U*D/SNU > -0 .457*ALDGi0(UST/W> 

F2=l.799-0.409*ALOG10(W*D/SNU) - 0 .314*ALOG10(UST/U) 
F3~ALOG10(ESP) 

E=Fi+F2*F3 
C=10.0**E 

CEE=C*S/(S-1. )*1 .0E-6 
1 CONTINUE 
RETURN 
END 

SUBROUTINE LAURSN(DFT,UST,S,V,GAMA,SUBW,SNU,DM,Y,G,W,CEE> 

THIS SUBROUTINE USES THE TOTAL-LOAD FORMULA DEVELOPED BY 
LAURSEN(19S8) 

D = DFT 
DP = DM 

CALL SHIELD<UST,DFT,SNU,SUBW,TC> 

UU=ALOG10(UST/W) 

IF(UU.GE.0.40) GO TO 1 

FF=1./(12.32-10.S*EXP(0.047*<UU+2.>)) 

GO TO 4 

IF<UU.GT.1.5) GO TO 2 
FF=2.04S*UU+0.942 
GO TO 4 

IF < UU,GE.2.2) GO TO 3 

FF=3.38+SQRT(1.416-(UU-2.51>**2 . ) 

GO TO 4 

FF=0.26*UU+3.953 
FL --1 0.0**FF 

TO=((GAMA*V*V)/(G*58.))*(DP/Y)**0.3333333 
IF(TO LT.TC) GO TO 5 

CEE = 0 01*( (TO/TO-l . )*FL*(D/Y)**1 . 1666666 
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GO TO 6 
CEE=0 . 0 
CONTINUE 
RETURN 
END 

SUBROUTINE SHIELD<UST,DFT,SNU,SUBU,TC> 

THIS SUBROUTINE COMPUTES THE CRITICAL BED SHEAR STRESS DERIVED 
FROM SHIELDS' FUNCTION 
REY=UST#DFT/SNU 
IF<REY.GT.10.0) GO TO 1 
TC = 0 .08*SUPW*DFT/REY**D.4 
GO TO 3 

IF(REY.GT.500.0 > GO TO 2 
TC = 0.022*SUBW*DFT*REY**0 . it 
GO TO 3 

TC=0.06*Sl)BU*DFT 

CONTINUE 

RETURN 

END 

SUBROUTINE SETVEL<D,T,S,W> 

THIS ROUTINE COMPUTES THE SETTLING VELOCITY OF SEDIMENT PARTICLES OF 
ANY DENSITY. THE ROUTINF. ASSUMES A SHAPE FACTOR OF 0 7 AND INTERPOLATES 
THE VALUES TABULATED BY THE SUBCOMMITTEE ON SEDIMENTATION, INTER¬ 
AGENCY COMMITTEE ON WATER RESOURCES, 1957. 

DIMENSION A (6,11),Z<2> 

DATA A<l,l),A<l,2),A(l,3>,A(l,4>,A(l,S>,A<i,6>,A<i,7>,A<l,8>, 

$A<1,9>,A(1,10),A<i,il)/ 

$0.04,0.06,010,0.20,0.40,0.80,1.50,200,3.00,700,10.00/ 

DATA A<2,1>,A<2,2),A<2,3>,A<2,4>,A<2,5>,A<2,6>,A<2,7),A<2,8>, 

♦A(2,9),A<2,10),A<2,11) / 

$0.10,0.24,0.60,1.80,4.60,9.50,16.10,19.90,25.30,39.50,44 00/ 

DATA A(3,1),A<3,2),A(3,3),A<3,4>,A(3,5),A(3,6),A(3,7>,A(3,8), 

$A < 3,9 >,A(3,10),A<3,11)/ 

$0.14,0.32,0 76,2.20,5.30,10.50,16.90,20.30,25.60,39 50,44.00/ 

DATA A<4,1),A<4,2>,A(4,3),A<4,4),A<4,S>,A(4,6),A<4,7),A(4,8), 

$A (4,9) , A < 4,1.0 > , A (4,11) / 

$0.18,0.40,0.92,2.50,5.80,11.00,17.SO,20.70,25.90,39 SO,44.00/ 

DATA A(5,1), A (5,2),A(5,3),A(S,4),A(5,S>,A<5,6),A<5,7),A<5,8), 

$A<5,9), A (5,10 ),A<5,11>/ 

$0.23,0.49,1.10,2.85,6.30,11.60,17.90,21.10,26.20,39.50,44.00/ 

DATA A<6,1),A<6,2), A < 6,3 ), A < 6,4 >,A<6,5>,A<6,6),A<6,7),A(6,8), 

$A<6,9),A(6,10),A<6,11>/ 

$0.29,0.57,1.26,3.20,6.70,12.0 0,18 10,21.50,26 SO,39.50,44.00/ 
VSC(T)=1.41E-5 -3.48E-7«(T-10.) +5.0OE-9*(T-l0.)*<T-15.) 

$ +2.67E-11*<T-10.>*<T —15 .>*(T-20.> -4 . OOE-12*(T-10.>* 

$ (T-15.)*<T-20.)*<T-25.> 

IF<D.GE.0.040) GO TO 2 
VISC=VSC(T) 

SS = 0. 0 009669*D*D/VISC 
GO TO 18 

INTERPOLATION 
CONTINUE 
Q=T/10. 
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KT=Q+i. 

PT=Q-KT+1. 

DL=ALOG10(D) 

DO 10 J=i,10 

IF(D.LE.A(1,J>) GO TO 20 
10 CONTINUE 

20 J=J-i 

C=ALQGi0(A < 1, J ) ) 

E=ALOG10<A<l,J+i>> 

PD=><DL-C)/<E-C> 

DO 50 L=i,2 
I=L+KT 

50 Z<L)=<i.-PD)*ALOGl0< A < I,J> >+PD*ALOGi0(A(I,J + i)> 
R=<1 -PT)*Z<i)+PT*Z<2> 

SS—10 HR 

C ADJUSTING SETTLING VELOCITY FOR SPECIFIC GRAVITY 

18 W=SS*SQRT<(S-l,0)/l.65)/30.5 

RETURN 
END 
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